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Rice '^‘University  has  developed  a  practical  computer  model  that  is 
capable  of  specifying  electron  and  ion  fluxes  in  the  middle 
magnetosphere  during  geomagnetic  storms.  The  model,  called  the 
Magnetospheric  Specification  Model  (MSM) ,  uses  ground-based  and 
satellite  data  from  the  Space  Forecast  Center-Environment  Data 
Base  to  establish  initial  and  boundary  conditions  and  to  determine 
(input  parameters  for  the  magnetic  and  electric  field  models.  These 
(input  values  are  updated  every  15  minutes,  and  new  output  fluxes 
jare  computed  for  the  same  times.  The  primary  function  of 
ithe  MSM  is  the  specification  of  fluxes  of  1-100  KeV  electrons  in 
the  geostationary  orbit  region.  However,  it  is  also  designed  to 
jspecify  a  broad  range  of  additional  parameters  for  the  global  iono- 
|spheric-magnetospheric  system,  including  fluxes  of  1-50  KeV  ions, 
jauroral  electron  precipitation  and  ionospheric  electric  fields. 
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The  model  is  accompanied  by  an  application  program  that  allows 
specification  of  fluxes  at  an  arbitrary  point  in  the  magnetosphere 
within  the  modeling  region.  Consistent  with  its  primary  function, 
the  MSM  has  been  tested  against  spacecraft  data  for  2  substantial 
storms  and  has  been  shown  to  produce  a  good  characterization  of  the 
enhancements  of  40  KeV  electron  fluxes  in  the  equatorial  plane. 
The  model  never  failed  to  predict  high  fluxes  when  they  were 
observed,  although  it  did  predict  high  fluxes  in  some  cases  when 
they  were  not  observed  and  it  did  fail  to  predict  flux  dropouts 
observed  by  the  spacecraft.  The  MSM  is  ready  for  adaptation  for 
usf»  in  an  operational  setting  wnere  the  goal  is  real-time  and 
retrospective  specification  of  hazardous  charged  particle  fluxes 
.associated  with  geomagnetic  storms. 
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1.  INTRODUCTION 
1.0.  Background 

The  magnetosphere  contains  plasma  and  energetic  particles  which,  during  times  of 
magnetic  disturbances,  pose  serious  hazards  for  spacecraft.  These  hazards  take  the  form  of 
spacecraft  surface  charging,  deep  dielectric  charging  and  disturbance  of  spacecraft 
electronics  by  direct  penetration  of  energetic  particles.  Most  commonly  the  effect  is  an 
electrical  transient  which  propagates  throughout  the  spacecraft  electronics.  Not  only  can 
temporary  or  permanent  malfunctions  result,  but  the  anomalies  may  mimic  intrinsic  or 
adversary  induced  electronic  component  failures. 

Diagnosis  and  prevention  of  environmentally  related  spacecraft  anomalies  requires  a 
knowledge  of  the  space  plasma  and  energetic  particle  fluxes  at  each  spacecraft  at  all  times. 
Since  most  satellites  do  not  carry  a  complete  suite  of  particle  measuring  instruments,  and 
even  if  they  did  they  could  not  forecast  local  conditions,  it  is  necessary  to  develop  other 
means  to  specify  the  state  of  space  weather.  Fortunately,  our  understanding  of  the  physics 
of  magnetospheric  processes  has  advanced  to  the  point  where  it  is  now  possible  to  build 
computer  models  which  can,  using  real-time  data  from  monitoring  spacecraft  and  ground 
observatories,  compute  particle  fluxes  at  arbitrary  positions  in  the  magnetosphere.  The 
Magnetospheric  Specification  Model  represents  the  first  such  operational  model  to  be 
developed.  As  this  contract  report  will  show,  this  effon  has  been  remarkably  successful. 

1.1.  Objectives 

The  objective  of  this  contract  has  been  to  develop  a  comprehensive  model  of  the 
Earth's  magnetosphere  for  operational  use  by  the  Air  Weather  Service.  The  model  is 
designed  to  specify  magnetospheric  conditions  using  real-time  data  provided  by  an 
environmental  database  at  the  Space  Forecast  Center.  The  computational  system  includes 
models  of  global  magnetic  and  electric  fields,  inner  magnetospheric  panicles,  and 
precipitating  auroral  particles.  It  can  be  operated  in  conjunction  wi  n  ionospheric, 
thermospheric,  solar-wind  and  solar-wind/magnetosphere  transpon  models  provided  by 
other  contractors.  The  model  makes  optimum  use  of  datasets  from  operational  Air  Force, 
DoE  and  NOAA  spacecraft  and  ground-based  observatories.  In  each  case,  a  back-up 
source  for  the  required  input  data  values  has  been  provided  in  order  that  the  model  can 
continue  to  function  in  the  absence  of  one  or  more  input  values. 
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2.  SCIENTIFIC  DESCRIPTION  OF  THE  MAGNETOSPHERIC 

SPECIFICATION  MODEL 


2.0.  Introduction 

This  scientific  description  is  intended  to  provide  an  overview  of  the  MSM.  More 
detail  on  certain  topics,  particularly  regarding  the  detailed  operation  of  the  MSM,  can  be 
found  in  other  sections  of  this  report  and  especially  in  the  Appendix. 

2.0.1.  Goals  and  Objectives 

The  goals  of  the  MSM  were  the  following: 

(i) .  It  should  produce  the  most  accurate  and  reliable  specifications  possible,  given  the 
available  input  parameters,  the  state  of  magnetospheric  science,  and  the  personnel  and 
computing  facilities  available  to  the  Air  Weather  Service. 

(ii) .  It  should  provide  a  limited  capability  for  forecasting,  specifically  upper-  and  lower- 
limit  forecasts. 

(iii) .  It  should  have  a  flexible  modular  structure,  to  allow  for  future  substitution  of  new 
and  improved  routines. 

(iv) .  It  should  be  comprehensive  enough  to  be  capable,  with  minor  alterations,  of 
specifying  a  wider  range  of  physical  parameters  than  is  now  required. 

(v) .  It  should  have  as  much  theoretical  consistency  as  possible,  given  the  constraints  listed 
under  item  1 . 

The  basic  objective  of  the  MSM  is  to  specify  certain  magnetospheric  parameters  based 
on  real-time  input  data  that  is  expected  to  be  available  to  the  Space  Forecast  Center.  The 
major  physical  parameters  calculated  by  the  MSM  include  the  following: 

(ii _ Magnetospheric  Particle  Fluxes. 

1-100  keV  Electrons  at  Synchronous  Orbit.  Calculation  of  these  fluxes  was  the 
primary  objective  for  the  MSM.  Test  data  were  made  available  to  us,  and  model 
predictions  have  been  extensively  checked  against  these  test  data,  as  discussed 
in  sections  3  and  4. 

Electrons  and  Ions  at  L  «  2  to  10.  Energies  100  eV  <  £  <  100  keV.  Calculation  of 
these  particle  fluxes  constituted  a  secondary  objective  for  the  MSM.  The  model 
calculates  these  fluxes,  but  very  little  data  have  been  available  to  allow  testing  of 
model  accuracy. 

Electrons  and  Ions  with  £  >  100  keV.  For  electrons  between  100  and  300  keV,  the 
MSM  estimates  fluxes  by  a  simple  algorithm  that  is  based  on  observed 
geosynchronous  fluxes.  For  electrons  above  3  MeV,  the  Koons-Gomey  model 
is  available  and  will  be  installed  shortly. 

(ii).  Energy  Fluxes  and  Average  Energies  of  Precipitating  Electrons  and  Ions.  These 
parameters  are  specified  by  the  MSM  and  are  thus  available  for  use  by  the  ionospheric 
model,  if  needed.  Comparison  of  predicted  electron  precipitation  against  data  for  April  21- 
23,  1988  is  discussed  in  Section  3.  The  statistical  electron-precipitation  model  of  [Hardy, 
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1985  #1].  Because  of  the  lack  of  real-time  data  for  testing  of  a  full  ion-precipitation 
algorithm,  the  MSM  is  set  to  use  the  statistical  ion-precipitation  algorithm  of  [Hardy,  1989 
#2]. 

(iiil.  Ionospheric  Electric  Fields.  The  ionospheric  potential  distribution  is  provided  for 
possible  use  by  the  ionosphere  model,  if  needed. 

(ivl.  Magnetospheric  E-  and  B-Fields.  These  fields  must  be  computed  by  the  MSM  as  part 
of  its  overall  modeling  of  the  magnetosphere.  The  B-field  model  is  directly  used  by  the 
application  program  MAPINT  to  map  fluxes  from  the  equatorial  plane  to  an  arbitrary  point 
on  a  field  line. 

2.0.2.  Logical  Structure  of  the  MSM. 

The  overall  logical  structure  of  the  MSM  is  diagrammed  in  Figure  2.1,  and  the  basic 
organizational  structure  of  the  presei.  section  is  based  on  the  figure.  Figure  2.2  illustrates 
the  major  subroutines  and  provides  a  more  detailed  flow  chart  of  the  MSM.  The  basic 
input  data  required  for  the  full  model  mn  are  read  from  the  Environmental  Database  and  are 
subjected  to  interpolation  operations  (described  in  Section  2.1.1).  In  most  cases,  the  input 
data  stream,  suitably  interpolated,  still  does  not  directly  specify  all  of  the  input  parameters 
required  by  the  MSM.  The  input  data  is  therefore  fed  into  a  set  of  front-end  models,  which 
utilizes  whatever  data  are  available  to  estimate  the  remaining  parameters  that  are  required  by 
the  central  routines  of  the  MSM  (Section  2.2).  Next  the  matrices  describing  the  essential 
characteristics  of  the  magnetic  field  are  computed  for  mark-time  during  the  period  to  be 
modeled.  (Fifteen-minute  intervals  have  been  used  in  our  tests.)  The  magnetic-field 
model,  which  represents  the  most  complex  element  of  the  MSM,  is  described  briefly  in 
Section  2.3,  and  fully  in  Appendix  f,  which  is  Robert  Hilmer's  Ph.D.  thesis.  After  the 
completion  of  the  computation  of  the  magnetic-field  matrices  for  the  modeled  period,  the 
program  computes  a  single  matrix  that  specifies  the  electrostatic  potential  at  every  grid 
point,  for  each  mark-time  during  the  event.  The  electric-field  model  is  described  in  Section 
2.4. 


2.1  Input  Data  from  the  Environmental  Database 

The  MSM  uses  real-time  data  from  a  variety  of  sources  to  establish  initial  and 
boundary  conditions  for  the  particle  distribution  functions,  and  to  establish  input 
parameters  for  the  E-  and  B-field  models.  These  data  are  retrieved  from  the  environmental 
database  through  a  subroutine  called  INDATA  (see  figure  2.2).  INDATA  is  called  by  the 
MSM  to  access  the  input  data  files  but  INDATA  must  be  written  by  the  contractor  preparing 
the  environmental  database.  (Note:  In  the  version  of  MSM  delivered  with  this  report, 
INDATA  exists  as  a  stub  to  access  the  test  datasets  delivered  with  the  model.)  As  input 
data  are  brought  in  they  are  tested  for  gaps  in  a  subroutine  called  PARGEN.  If  data  are 
missing  from  a  given  input  parameter  for  a  time  gap  in  excess  of  a  certain  time,  "front-end" 
models  are  used  to  generate  substitute  values  for  the  missing  data  values.  The  front-end 
models  rely  on  empirical  algorithms  described  in  section  2.1.3.  The  data  are  also 
interpolated  to  obtain  values  at  the  model  step  times.  This  interpolation  is  described  in 
section  2.1.2. 

In  addition  to  the  data  brought  in  through  INDATA,  that  are  used  as  input  to  the 
models  primary  calculations,  the  MSM  also  uses  geostationary  satellite  particle  flux  data  to 
provide  a  comparison  with  the  model  output. 
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2.1.0.  MSM  Input  Parameters 

The  MSM  requires  seven  basic  input  parameters,  plus  two  sets  of  geostationary 
particle  data  for  estimates  of  errors  in  the  model  ouqjut.  These  parameters,  together  with 
their  applications  in  the  MSM  their  variable  names  and  units  are  as  follows: 


Input  Parameters  Used  by  the  MSM 

Table  2.1 


Input  parameter  and  use 

Name 

Unit 

Kp  (Used  for  initial  and  boundary  fluxes) 

FKP 

NCTE 

Dst  (Determines  ring  current  strength  in  B-field  model) 

DST 

NANOTESLA 

LOW  LAT  BOUNDARY  OF  THE  AURORA  (Determines  auroral 
zone  boundary  in  E-field  model  and  is  used  to  constrain  the 
B-field  mapping) 

DLATAZ 

DEGFES 

POLAR  CAP  POTENTIAL  (Input  to  E-field  model) 

PCP 

KILOVOLTS 

POLAR  CAP  POTENTIAL  PATTERN  (Input  to  E-field  model) 

XIPATT 

NCfE 

SOLAR  WIND  VELOCITY  (Determines  standoff  distance  in  B- 
field  model) 

SWVEL 

Km/g 

SOLAR  WIND  DENSITY  (Determines  standoff  distance  in  B- 
field  model) 

SWDEN 

PROT/c^3 

GEO  ELECTRON  FLUX  (Used  to  check  model  predictions)  * 

EPSAT1 

^‘■^^^CM^-S-Sr-Kev 

GEO  ION  FLUX  (Used  to  check  model  predictions)  * 

EPIONS 

'°^®/cM2.S-S,Kev 

*  -  Not  loaded  by  calls  to  INDATA. 

The  low  latitude  boundary  of  the  aurora  is  determined  from  the  electron  precipitation 
data  from  the  DMSP  satellites  through  an  algorithm  developed  by  the  Geophysics 
Laboratory.  The  boundary  position  is  adjusted  for  midnight  magnetic  local  time  using  the 
algorithm  published  by  Gussenhoven  et  al.[1983]. 

The  polar  cap  potential  drop  is  derived  from  DMSP  ion  drift  meter  data  by  an 
algorithm  developed  at  the  University  of  Texas  at  Dallas  by  Rod  Heelis  and  Marc  Hairston. 

The  polar  cap  potential  pattern  type  is  from  Maynard  -  Heppner  as  digitized  by  Fred 
Rich.  Three  patterns  are  possible.  The  choice  of  type  is  determined  by  an  algorithm 
provided  by  the  University  of  Texas  at  Dallas  group.  For  the  tests  run  to  date,  only  two 
pattern  types  have  been  required.  These  are  representative  of  the  IMF  conditions  Bz  south 
and  north. 

Solar  wind  data  is  intended  to  be  from  a  solar  wind  monitor.  For  our  tests  this  data 
was  derived  from  the  NSSDC  OMNI  dataset  (King,  1989], 

The  GEO  electron  and  ion  fluxes  are  obtained  from  geostationary  satellites  whenever 
possible  and  are  read  by  the  MSM  differently  from  the  rest  of  the  parameters. 


Kp  and  Dst  are  provided  by  the  Air  Force  and  NOAA. 
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2. 1.0.1.  Parameters  Called  by  INDATA 

Each  input  parameter  (except  the  GEO  electron  and  ion  fluxes)  is  loaded  separately  by 
calling  INDATA.  .  The  form  of  the  call  to  INDATA  is:  CALL  INDATA(param- 
name,STARTT^NDT,NDIMT)ARRY,NUMNUM) 

The  argiunents  of  the  call  have  the  following  meanings: 
param-name  is  a  valid  parameter  name  from  the  MSM  input  parameter  list. 

STARTT  is  the  starting  time  for  the  current  data  request. 

ENDT  is  the  ending  time  for  the  current  data  request 

NDIM  is  the  maximum  number  of  input  parameter  values  that  may  be  returned  and  the 
horizontal  dimension  of  the  data  array,  DARRY. 

PARRY  is  the  two  dimensional  array  to  be  retrieved.from  the  environmental  database.  The 
size  of  DARRY  will  be  NDIM  by  8. 

NUMNUM  is  the  number  of  data  values  placed  in  DARRY  by  INDATA. 

The  format  consists  of  an  eight-word  record  as  follows: 

Word  1  =  data  value  at  t 

Word  2  =  year  at  t  (full  year  required  eg.  1990) 

Word  3  =  decimal  day  at  t 

Word  4  =  spacecraft  geomagnetic  latitude  at  t 

Word  5  =  spacecraft  geomagnetic  longitude  at  t 

Word  6  =  spacecraft  altitude  at  t 

Word  7  =  magnetic  local  time  at  t 

Word  8  =  data  error  quality  value  at  t 

where  t  represents  the  time  tag  for  the  data  value,  ie.  words  2  and  3  and  data  values  are 
expected  one  after  the  other  in  time-sequential  order  with  a  complete  block  of  all  available 
data  making  up  one  call  to  INDATA. 

The  last  five  words  are  ignored  in  the  present  version  of  the  MSM. 

Data  words  not  used  or  empty  are  represented  by  -999. 

All  parameters  are  assumed  to  be  single  precision,  floating  point  numbers. 
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2.1.02.  The  GEO  Satellite  Particle  Data 

The  electron  and  ion  flux  data  from  geostationary  satellites  are  not  used  to  set  model 
conditions  and  Uiey  must  be  formated  with  greater  flexibility  because  of  the  variability  of 
detector  energy  channels  available  and  the  possibility  of  further  changes  in  detector  suites  in 
the  future.  Therefore,  these  data  are  not  handled  in  the  same  manner  as  the  other  seven 
input  data  parameters.  They  are  instead  called  directly  as  files  with  the  the  following 
format: 


Header  record: 


Number  of 

Year 

Decimal  day 

Geographic 

Geographic 

Altitude 

energy 

(eg. 

(eg.  112.500 

latitude  of 

longitude  of 

of 

local  time 

channels  in 

1988) 

=noon,  April 

satellite 

satellite 

satellite 

of  the 

the  following 
data  record 

21  in  leap  yr) 

(degrees) 

(degrees- 

east) 

(Km) 

satellite 

(HrMin.) 

Data  record: 


Low  end  of 
energy  channel 

1  (KeV) 

High  end  of 
energy  channel 

1  (KeV) 

Log  10  Rux 
(No./cm2-S-Sr- 
KeV)  channel  1 

1  sigma  error  in 
Flux  channel  1 

Species 

(-l=electrons, 

l=ions) 

Low  end  of 
energy  channel 
2  (KeV) 

High  end  of 
energy  channel 
2  (KeV) 

Log  10  Flux 
(No./cm2-S-Sr- 
KeV)  channel  2 

1  sigma  error  in 
Flux  channel  2 

Species 

(-l=electrons, 

l=ions) 

Low  end  of 
energy  channel 
n  (KeV) 

High  end  of 
energy  channel 
n  (KeV) 

Log  10  Flux 
(No./cm2-S-Sr- 
KeV)  channel  n 

1  sigma  error  in 
Flux  channel  n 

Species 

(-l=electrons, 

l=ions) 

All  variable  are  assumed  to  be  single  precision  floating  point  numbers. 

Data  words  not  used  or  empty  are  represented  by  -999. 

Data  from  different  satellites  may  be  interleaved  but  are  assumed  to  be  time 
sequential. 

2.1.1.  Input  Parameter  Interpolation 

Input  data  values  have  time  tags  that  not  generally  simultaneous  with  the  time  steps 
for  a  given  model  run.  (The  model  nominal  time  step  is  a  1 5  minute  interval  beginning 
with  the  model  start  time.)  Since  input  data  values  are  required  at  each  model  time  step,  a 
linear  interpolation  is  performed  between  nearest  neighbor  data  values.  Several  methods  of 
smoothing  the  data  were  examined  and  rejected  because  they  resulted  in  a  loss  of  time 
resolution  or,  in  the  case  of  spline  fits,  caused  unrealistic  overshoot  (see  Quarterly  Status 
Report  no.  10.  The  required  interpolation  occurs  in  the  subroutines  DNTTRP  and 
DTXIPT. 

At  same  time  that  the  interpolation  takes  place,  a  time  gap  sensing  routine  determines 
data  gaps  between  the  required  model  time  step  time  and  the  nearest  available  data  time. 
This  is  done  so  that  significant  data  gaps  can  be  filled  by  proxy  values  for  the  input 
parameters.  These  proxy  "values"  are  determined  by  the  Front-end  models.  The  maximum 
allowable  time  gap  before  proxy  values  are  used  was  determined  separately  for  each 
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parameter  by  a  statistical  study.  In  this  study,  for  a  run  of  data  from  April  21-23,1988,  an 
increasing  number  of  data  points  were  intentionally  deleted  and  replaced  by  linear 
interpolations  and  the  RMS  error  computed  after  each  new  deletion.  Ifroxy  data  values 
were  then  inserted  in  place  of  the  interpolations  and  the  RMS  error  computed  again.  At 
each  point,  the  two  RMS  errors  were  compared.  When  the  RMS  error  obtained  with  proxy 
values  inserted  became  lower  than  that  for  no  data  values  but  with  interpolations,  the 
maximum  time  gap  was  considered  to  be  established. 

2.2.  Front>end  Models  and  Default  Models 

An  important  constraint  placed  on  the  MSM,  early  in  the  program,  was  that  it  must 
"fail  gracefully".  This  meant  that  the  model  must  maintain  it's  accuracy  and  validity  in  the 
event  of  loss  of  some  or  even  all  of  the  input  data  parameters,  except  Kp.  In  order  to  meet 
this  constraint,  two  strategies  were  devised. 

The  first  strategy  was  the  development  of  Front-end  models  which  provide  values  for 
the  input  data  parameters  if  one  of  the  values  were  missing  from  the  input  data  stream, 
either  intermittendy  or  for  a  long  period.  A  separate  Front-end  model  is  required  for  each 
input  data  parameter  and  in  some  cases  several  backup  models  are  required. 

The  second  strategy  involved  the  use  of  so-called  "Default"  models  to  completely 
replace  the  MSM  computed  flux  output  if  the  only  remaining  input  data  parameter  is  Kp. 

2.2.1.  Front-end  Models 

As  indicated  above,  data  gaps  in  the  real-time  input  parameters  must  be  filled  with 
proxy  values.  Front-end  models  which  rely  on  other  parameters  present  in  the  data  stream 
were  developed  to  provide  these  values.  In  one  case,  the  parameter  required  by  the  model 
itself  is  a  secondary  parameter  derived  from  the  input  parameters,  solar  wind  velocity  and 
solar  wind  density.  This  parameter  is  the  standoff  distance. 

Since  it  was  anticipated  that  the  Front-end  models  would  be  primarily  empirically 
based  and  would  probably  be  dependent  on  Kp.  Dst,  or  similar  geophysical  parameters,  the 
first  step  was  an  extensive  literature  search  to  look  for  paper  relating  various  parameters. 
This  search  was  conducted  by  Robert  Hilmer  and  took  a  one  month.  The  search  netted  131 
papers.  These  are  listed  in  the  Appendix  as  the  MSM  Reference  List:  Correlation  Studies. 
The  papers  were  arranged  in  a  matrix  format  such  that  all  papers  relating  any  two 
parameters  could  be  quickly  identified.  This  matrix  is  also  reproduced  in  the  Appendix. 

Next,  algorithms  were  developed  which  made  use  of  the  relationships  found  in  the 
literature.  In  one  case,  published  relationships  were  put  together  to  yield  a  new  relationship 
between  the  low  latitude  boundary  of  the  aurora  and  the  inner  edge  of  the  plasma  sheet.  In 
another  case,  a  new  empirical  relationship  was  established  between  Kp  and  Dst. 
Development  of  these  Front-end  model  algorithms  is  discussed  in  detail  in  Quarterly  Status 
Reports  nos.3,  9,and  10. 

The  following  table  indicates  the  sources  of  the  input  parameter  values  and  references 
to  addition  input  values  required  by  the  MSM: 


Table  2.2 
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FRONT  END  MODEL  INPUT  PARAMETER  SOURCES 


PARAMETER 

USED  IN 

METHOD  OF  DETERMINATION  (in  order  of  priority) 

Dst 

B-field 

1 .  Ground  magnetometers 

2.  Dst=  -  355  +5.25  *  DLATAZ 

3.  Dst=  -Kp*4  for  Kp<4  :  Dst=  -20-((Kp*10)-40)*2.6  for 
Kp>4 

4.  Persistence  if  less  than  6  hours 

Low  latitude  B-field  1 .  DMSP  or  NOAA  satellite  data 

boundary  of  the  2.  DLATAZ=66.95  -  2.03  *  Kp  (Gussenhoven  et  al.,  1983) 

aurora 


Polar-cap 
potential  drop 


Polar-cap 
pattern  type 


E-field  1.  Ion  drift  data  from  DMSP  (Heelis  and  Hairston). 

2.  PCP(Kv)=1 4.587  +  1.7  *  (Kp*10)  (Reiff,  Private  comm.) 


E-field  1.  Ion  drift  data  from  DMSP  (Heelis  and  Hairston) 
2.  Default:  assume  pattern  #1 


Standoff  B-field  1.  If  solar  wind  parameters  are  available,  then: 

Dist.(replaces  rs  =  f(f  Bo2)/2  %  po  where  f=1.35,  Bo=surf.  field  at 

SWVEL&  the  equator,  po*S.W.  mass  density,  and  Vo=  the  S.W.  velocity, 

SWDEN  if  otherwise: 

2.  Rs=  I1.y-Kp,  but,  Rs  >  or  =  6  Rp  at  ail  times 


Inner  edge  of  B-field  See  section  2.3 
plasma  sheet  at 
midnight  in 
eauatorial  plane 


Poleward  edge  E-field  See  section  2.4.2 
of  model  region 


Poleward  edge 
of  field- 

E-field 

See  section  2.4.2 

reversal  region 

Partrac  See  section  2.5 


Initial  p,Ti,Te  Partrac  See  section  2.5 


Notice  that  the  lowest  priority  algorithm  is  either  Kp  based  or  else  a  fixed  default 
value.  This  means  that  the  MSM  can  run  on  Kp  alone. 
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2.2.2.  Default  Models 

In  order  to  meet  the  constraint  that  the  MSM  produce  useful  output  under  conditions 
when  only  Kp  is  available  as  input  data,  two  Kp  based  statistical  models  were  incorporated 
to  replace  the  computed  output  of  the  MSM  in  the  event  of  loss  of  all  input  data  except  Kp. 
These  were  the  Hardy  et  al.,[1985]  model  of  precipitating  electron  and  ion  fluxes  in  the 
auroral  zone  and  a  model  developed  by  Garrett  for  the  electrons  and  ions  at  geostationary 
orbit. 


It  was  determined,  during  testing  that  the  MSM  run  on  Kp  alone  produces  better 
output  than  the  Garrett  model,  however,  the  Garrett  model  has  proven  useful  because  it  is 
used  to  establish  the  initial  conditions  for  the  model  in  the  vicinity  of  the  geostationary 
orbit.  It  is  probable  that  the  MSM  computed  auroral  precipitation  output  based  on  Kp  is 
also  better  than  the  Hardy  et  al.  model  output,  but  we  have  not  been  able  to  confirm  this 
since  DMSP  precipitating  particle  data  had  not  been  made  available  by  the  end  of  the 
contract.  Both  the  Hardy  et  al.and  Garrett  models  are  average  models  and  therefore  contain 
no  intrinsic  time  dependence  other  than  that  derived  from  the  variation  in  Kp.  We  also 
evaluated  for  possible  application,  an  average  auroral  precipitation  model  by  Foster  et 
al.[1986].  This  model  was  not  used  because  it  is  driven  by  a  precipitation  power  index  not 
re^ly  available. 

2 2 2.1.  The  Hardy  et  al.  Model  of  Precipitating  Auroral  Fluxes 

The  Hardy  et  al.model  is  actually  two  models,  one  for  electrons  and  one  for  ions. 
They  are  represented  in  the  MSM  by  the  subroutines  AURLl  and  AURL2S  for  electrons 
and  ions  respectively. 

The  electron  model  was  published  by  hardy  et  al.  [1985]  and  updated  coefficients 
were  provided  to  us  by  W.  J.  McNeil  of  Radex  Inc.  and  David  Hardy  of  GL.  The  model 
employs  an  Epstein  function  to  model  the  auroral  electron  fluxes  in  latitude  and  a  Fourier 
series  for  the  local  time  dependence.  The  Epstein  function  has  the  form 


e(h)=r+s,(h-hj+(s2-s,)  x  In 


(2.2-1) 


where  e(h)  is  the  energy  or  number  flux  dependent  on  the  peak  flux,  r,  the  latitude,  h,  the 
latitude  of  peak  flux,  ho,  and  si  and  S2,  the  slopes  on  either  side  of  the  peak.  The  Fourier 
series  has  the  form 


a(T)=2^  CnCos{-jy-)-i-S„sin(-y^) 

(2.2-2) 

where  a  is  the  calculated  parameter  which  can  be  either  r,  ho,  si,  or  S2.  For  each  Kp 
interval,  this  analytic  function  requires  52  coefficients  represented  in  the  above  equation  by 
C  and  S.  These  coefficients  are  embedded  in  the  MSM  and  can  be  found  in  the  appendix. 
The  model  is  used  as  follows:  to  determine  an  energy  flux  for  a  given  latitude,  h,  and  local 
time,  T,  use  the  Fourier  series  to  determine  r,  ho,  si,  and  S2.  Then  insert  these  into  the 
Epstein  function  along  with  h. 
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22.22  The  Garrett  Model  of  Geostationary  Orbit  Fluxes 

Garrett  et  al.  [1981a,  1981b]  published  a  statistical  analysis  of  the  geostationary  orbit 
electron  and  ion  environment  based  on  data  from  ATS  5  and  6.  Garrett  and  DeForest 
[1979]  had  previously  outlined  an  analytical  approach  to  modeling  such  average  flux  values 
in  terms  of  the  moments  of  the  distribution  functions.  Early  in  the  MSM  project  it  was 
realized  by  John  Gaudet  that  this  early  work  formed  the  nucleus  upon  which  could  be  built 
a  more  substantial  characterization  of  the  geostationary  fluxes  which  would  be  useful  to  the 
MSM.  Based  on  this,  using  additional  data  from  Scatha  and  P-78,  in  addition  to  the  ATS  5 
and  6  data.  Hank  Garrett  and  John  Gaudet  prepared  a  new  model  (or  rather  set  of  models) 
of  geostationary  electrons  and  ions  for  the  MSM.  Akira  Nagai  then  adapted  these  for  use  in 
the  MSM.  Generically  these  are  refered  to  as  the  Garrett  model  and  they  are  embedded  at 
numerous  places  in  the  MSM  subroutines. 

The  Garrett  model  uses  a  Fourier  expansion  of  the  diurnal  and  semi-diurnal  local  time 
components  and  a  linear  fit  to  Ap.  The  fits  determine  the  coefficients  of  the  expansions  of 
the  four  moments  of  the  distributions  functions.  These  moments  are  then  convened  to  bi- 
maxwellian  distribution  functions  which  are  integrated  to  provide  the  fluxes  in  any  required 
energy  passband.  The  equation  for  the  moments.  Mi  is 


Mi=  Cl C2COS  (^)+C3sin  (i^) 
C,cos(l^)+C5sin(i^)-h 
Ap  [C6+  C7COS  (^)+C8sin  (1^)  + 
C9COS  (^^)+Ciosin  (^^)] 


(2.2-3) 


The  Cs  are  the  coefficients  obtained  by  Garrett  and  Gaudet,  and  It  is  the  local  time.  Ap  is  8 
times  the  linearized  value  of  Kp  times  3.  There  are  eight  moments,  M),  number  density  and 
number  flux  and  energy  density  and  energy  flux  for  the  electrons  and  ions. 

During  testing  of  the  Garrett  model,  it  was  determined  that  the  there  were  problems 
with  the  model  predicting  unphysically  high  temperatures  at  high  Kp  near  midnight.  This 
was  solved  by  restricting  Kp  to  values  of  6  or  lower  when  the  Garrett  Model  is  used. 
Also,  the  Garrett  model  gives  unrealistically  low  fluxes  at  higher  energies  because  of  the 
low  high-energy  tail  inherent  in  the  maxwellian  distribution. 
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2.3.  The  Magnetic-Field  Model 
2.3.0.  Introductory  Comments 

Over  the  course  of  contract  F19628-87-K-0001,  a  completely  new  magnetic-field 
model  was  developed  for  use  in  the  Magnetospheric  Specification  Model  and  the  Rice 
Convection  Model.  This  new  magnetic-field  model,  which  constitutes  the  most  complex 
component  of  the  MSM,  represents  several  years  of  work  by  Robert  Hilmer,  carried  out 
under  the  supervision  of  Hannes  Voigt.  We  attempt  only  a  brief  summary  of  this  model 
within  the  main  body  of  this  repon.  A  complete  and  detailed  description  is  contained  in 
Hilmer's  Ph.D.  thesis,  which  is  included  as  Appendix  f 

The  major  goal  in  this  endeavor  was  a  B-field  model  that  would  be  flexible  enough  to 
represent  a  wide  variety  of  magnetospheric  B-field  configurations,  using  input  parameters 
that  would  be  available  in  the  MSM  input  datastream.  The  observation-based  B-field 
models  that  are  most  commonly  used  by  magnetospheric  physicists  today  are  those  of 
Tsyganenko  [1987,1989]  and  Tsyganenko  and  Usmanov  [1982],  but  those  models,  which 
are  solidly  based  on  averaged  magnetic-field  observations,  characterize  the  state  of  the 
magnetosphere  by  a  single  parameter  Kp.  However,  the  magnetic-field  configuration  can 
vary  widely  for  a  given  vdue  of  Kp.  The  MSM  is  designed  to  run  from  much  more 
detailed  red-time  observational  information  than  just  Kp.  Therefore,  a  magnetic-field 
model  was  needed  that  could  respond  flexibly  to  independent  variations  in  several  input 
parameters. 

2.3.1.  Inputs  to  the  Magnetic-Field  Model. 

The  input  parameters  used  by  the  MSM's  ma^etic-field  model  are  the  following: 

(i)  Magnetopause  standoff  distance,  which  is  estimated  from  the  solar-wind  pV^  if  solar- 
wind  data  are  available.  Otherwise,  the  standoff  distance  is  estimated  from  Kp  and  a 
statistical  relation.  (See  Section  2.2.). 

(ii)  Qjst,  which  is  derived  from  ground  magnetometers  stations. 

(iii)  Mapping  of  the  equatorward  edge  of  the  diffuse  aurora  at  local  midnight.  The 
latitude  of  the  equatorward  edge  is  estimated  from  real-time  DMSP  electron  data,  using  an 
algorithm  supplied  separately  by  the  USAF  Geophysics  Lab.  We  have  used  this  direct 
observational  information,  along  with  published  statistical  information,  to  place  a  restriction 
on  the  mapping  of  the  field  line  that  connects  the  equatorward  edge  of  the  aurora  at  local 
midnight  to  the  inner  edge  of  the  plasma  sheet.  The  description  presented  here  supersedes 
that  presented  in  Section  5.4  of  Appendix  f.  We  know  of  no  study  based  on  observational 
statistics  that  relates  h.mid^  the  invariant  latitude  of  the  equatorward  edge  of  the  diffuse 
aurora  at  local  midnight,  to  r^,  the  equatorial  geocentric  distance  of  the  inner  edge  of  the 
plasma  sheet  at  the  same  local  time.  However,  statistical  studies  have  related  each 
parameter  to  Kp.  The  results  of  Gussenhoven  et  al.  [1983]  suggest  that 

^mid  =  66.95°  -  2.03  ATp  (2.3-1) 

Statistical  analysis  of  the  /STp-dependence  of  r^id  is  somewhat  more  difficult.  However, 
we  have  used  the  published  statistical  studies  of  Horwitz  et  al.  [1986]  and  Kivelson 
[1976]  to  obtain  rough  upper- and  lower-estimate  linear  formulas  relating  r^id,  to  Kp  .  To 
get  an  idea  of  the  relationship  between  Ao  and  r^id  we  eliminated  Kp  between  those  two 
formulas  ^d  equation  (2.3.1).  The  resulting  upper-  and  lower-estimate  formulas  for 
as  a  function  of  Amid  are  displayed  as  rhigh  and  rlow  in  Figure  2.3.1.  The  data  on  which 
those  formulas  are  based  came  mainly  from  the  geosynchronous-orbit  region  of  the 
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magnetosphere,  and  we  expect  those  formulas  to  lose  validity  when  the  equatorward  edge 
is  outside  the  range  from  4  to  7  /?£.  It  is  clear,  in  fact,  that  the  formulas  give  ridiculous 
values  for  Amid  ^  48°.  However,  it  is  useful  to  bear  in  mind  another  simple  approximate 
relationship  between  rmid^  and  A^,  namely  the  dipole  formula: 

^ dipole  ~  [cos(A^id)]-^  (2.32) 

which  is  also  displayed  in  the  figure.  We  expect  that  r^id  will  always  exceed  rdipou^ 
because  the  midnight  field  is  normally  expected  to  be  inflat^  relative  to  the  dipole,  but  we 
also  suspect  that  the  ratio  r„ud  I r dipole  will  approach  unity  when  the  inner  edge  comes  close 
to  the  Earth  in  times  of  extremely  high  activity,  because  precipitation  and  charge-exchange 
loss  near  the  Earth  should  rapidly  reduce  the  convected  plasma-sheet  population  and  inhibit 
its  capability  to  distend  the  field  lines.  Another  extreme  case  with  which  the  MSM  will 
occasionally  have  to  deal  is  that  in  which  the  measured  equatorward  edge  of  the  aurora  is 
indicated  to  be  at  extremely  high  latitude.  It  is  our  understanding  that  the  inner  edge  of  the 
plasma  sheet  at  local  midnight  is  rarely,  if  ever,  observed  to  be  at  >  10  Re  geocentric 
distance.  With  all  these  considerations  in  mind,  we  chose  for  the  MSM  the  following 
prescription  for  the  mapping  of  the  low-latitude  edge  of  the  aurora  at  local  midnight: 


^mid  _  1  + 

dipole^ 

for  r dipole  ^  5 

(2.3.3fl) 

—  It 

' dipole 

1  5  ) 

=  125 

l' dipole 

for 

5  ^  r dipole  ^  8 

(2.3  3&) 

O 

II 

3 

for 

^ dipole  ^  8 

(233c) 

Equation  2.3.3  does  not  appear  in  the  main  MSM  or  the  MAPINT  program  anywhere.  It 
was,  however,  used  continually  in  the  generation  of  the  magnetic-field  matrices  at  Rice.  It 
is  described  here,  since  it  represents  a  significant  approximation  that  should  be  borne  in 
mind  when  the  MSM  is  used. 

fivl  (Collap.sel.  The  Hilmer-Voigt  model  includes  a  feature  that  allows  collapse  of  the 
midnight  region  of  the  magnetotail,  as  apparently  occurs  in  the  expansion  phase  of  a 
magnetospheric  substorm.  Although  this  is  an  interesting  feature  from  a  scientific 
viewpoint,  we  are  not  utilizing  it  in  the  present  version  of  the  MSM,  because  we  have  not 
found  an  automatic  algorithm  for  deducing  substorm-associated  tail-field  collapse  in  the 
data  stream  that  is  expected  to  be  routinely  available  for  MSM  operations. 
fvl  (Dipole  tilt).  The  Hilmer-Voigt  model  allows  the  Earth's  magnetic  dipole  to  be  tilted 
away  from  the  nominal  direction  perpendicular  to  the  solar-wind  velocity  vector. 
However,  the  way  in  which  the  magnetic-field  model  has  had  to  be  incorporated  within  the 
MSM  —  in  terms  of  pre-computed  matrices  --  has  precluded  using  this  basic  feature  of  the 
B-field  model. 

2.3.2.  Pre-Computed  Magnetic-Field  Matrices. 

The  computing  facilities  that  are  now  expected  to  be  available  for  routine  operation  of 
the  MSM  preclude  real-time  computation  of  magnetic-field  models.  Therefore,  we  have 
precomputed  a  super-grid  of  magnetic-field  models,  in  which  each  model  is  characterized 
by  the  following  three  parameters: 

(i)  Magnetopause  standoff  distance  (6,  8, 10,  12,  and  14  Re). 

(u)  Dst  (-400,  -300,  -200,  -150,  -100,  -50,  0,  and  50  nT). 
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(iii)  Equatorward  edge  of  the  diffuse  aurora  at  local  midnight  (16  latitudes  from  49.5°  to 
69.3°). 

Each  B-field  model  in  this  supergrid  is  characterized  by  five  matrices:  XMIN,  YMIN, and 
ZMIN,  which  specify  the  location  of  the  "equatorial  mapping  point"  of  the  field  line  in 
GSM  coordinates,  for  each  point  in  the  fixed  ionospheric  grid;  BMIN,  the  equatorial  field 
strength;  and 


VM  =  [!ds/B]-^^^  (2.3.4) 

where  the  integral  extends  along  a  field  line,  from  the  southern  ionosphere  to  the  northern. 
The  parameter  VM  determines  the  gradient/curvature  drift  of  the  particles.  The  five 
matrices  provide  all  necessary  magnetic-field  information  for  running  the  main  MSM.  We 
interpolate  on  the  supergrid  of  matrices  to  find  the  five  matrices  that  represent  the  essential 
characteristics  of  the  magnetic  field  at  any  given  time. 

It  should  be  remarked  that  two  of  the  above-mentioned  matrices,  ZMIN  and  BMIN, 
are  not  actually  used  in  the  present  version  of  the  MSM.  They  are  being  included  to 
facilitate  future  enhancement  of  model  capabilities:  ZMIN  will  be  required  if  dipole  tilt  is 
included  in  future  versions  of  the  model;  BMIN  is  included  to  facilitate  explicit 
computation  of  drift  velocities  for  particles  that  mirror  near  the  equatorial  plane,  should  that 
capability  be  needed  in  special  cases. 

An  application  program,  MAPINT,  can  map  from  an  arbitrary  point  P  within  the 
modeled  region  of  the  magnetosphere  to  the  equatorial  plane  or  ionosphere,  and  thus  locate 
point  P  within  our  grid.  This  application  program  is  described  in  Section  3. 

The  Hilmer-Voigt  magnetic  field  is  computed  as  a  superposition  of  the  fields  due  to 
five  basic  currents,  as  follows: 

(i)  The  Earth's  main  field,  currently  assumed  to  be  a  dipole  B^; 

(ii)  the  field  due  to  the  ring  current; 

(iii)  B,au,  the  field  due  to  the  cross-tail  current; 

(iv)  the  field  due  to  the  Chapman-Ferraro  currents,  which  flow  of  the  magnetopause. 

The  technical  details  of  the  Hilmer-Voigt  magnetic-field  model  are  described  in 
Appendix  f. 
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2.4.  The  Eiectric-Field  Model 
2.4.0.  Introductory  Comments. 

A  new  model  of  the  ionospheric  electric  field  was  developed  for  the  MSM,  a  model 
that  was  specifically  designed  to  be  driven  by  observational  parameters  that  will  be 
available  in  real-time  for  the  MSM,  specifically  the  polar-cap  potential  drop,  the  polar-cap 
pattern  type,  and  the  equatorward  edge  of  the  diffuse  aurora  at  local  midnight. 

2.4.1.  Division  of  the  Ionosphere  into  Regions. 

The  essence  of  the  MSM's  electric-field  model  involves  dividing  the  ionosphere  up 
into  four  different  regions  as  shown  in  Figure  2.4-1.  The  four  regions,  which  are  treated 
quite  differently,  are  the  following: 

Region  0:  The  polar  cap,  inside  the  boundary  marked  a  above,  which  is  the  poleward 
boundary  of  the  "electric-field-reversal  region"; 

Region  1:  The  electric-field-reversal  region;  within  our  simplified  view,  this  is  coincident 
with  the  area  occupied  by  the  region- 1  Birkeland  currents; 

Region  2:  The  region  of  the  main  auroral  sunward  flow;  it  extends  from  boundary  b, 
which  is  the  equatorward  edge  of  the  electric-field  reversal  region,  to  boundary  c,  which  is 
the  equatorward  edge  of  the  shielding  layer; 

Region  3:  The  low-latitude  region,  from  c  to  the  magnetic  equator. 

The  three  boundaries  (a,  b,  and  c)  are  defined  to  be  ellipses, 
following  equation: 

(x-DXf  ^  (y-DYf  ^  j 

where 

X  =  Q  cos((|)) 
y  =  9  sin((I)) 

where  0  is  invariant  colatitude,  and  ())  is  magnetic  local  time  (0  at  noon,  n/2  at  dusk...). 
We  do  not  make  any  a  priori  assumptions  about  whether  A  or  B  is  bigger. 

2.4.2.  Determination  of  the  Ellipse  Parameters  A,  B,  DX,  and  DY. 

In  the  program,  the  parameters  defining  the  three  electric-field  ellipses  shown  in 
Figure  2.4-1  are  defined  by  the  three-component  vectors  4(3),  B(3),  DXO),  and  DYO). 
In  each  case,  the  first  component  represents  boundary  a,  the  polewardmost  boundary  that 
lies  at  the  high-latitude  edge  of  the  electric-field-reversal  region,  component  2  represents 
boundary  b,  which  lies  at  the  equatorward  edge  of  the  electric-field-reversal  region,  and 
component  3  represents  boundary  c,  which  lies  at  the  equatorward  edge  of  the  shielding- 
layer  region. 

2. 4. 2.1.  Location  of  Ellipse  2,  i.e.,  Boundary  b. 

We  have  chosen  to  associate  this  boundary,  which  in  the  ionosphere  is  the 
equatorward  edge  of  the  electric-field-reversal  region,  with  the  high-L  boundary  of  our 
main  modeling  region.  Furthermore,  we  choose  the  following  "reasonable"  locations  for 
this  boundary  in  the  equatorial  plane: 


which  satisfy  the 
(2. 4- la) 

(2.4-1/?) 

(2.4-lc) 


Section  2,  Page  14 


r  moddl  bcufid.  — 

^ noon  “ 

0.95  r standoff 

i2.4-2a) 

rdawn'^'^"^'^  =  standoff 

(,2.4-2b) 

r  •  j  •  L  .ffutdil  bchuuL 
^  midnight 

—  2.00  r standoff 

(2.4-2C) 

Specifically,  we  define  bound,  jq  5^  thg  equatorial  geocentric  distance  of  a  field 

line  that  lies  at  the  equatorward  edge  of  the  electric-field-reversal  region  and  crosses  the 
conducting  layer  of  the  ionosphere  at  local  dawn.  It  will  not  cross  the  equatorial  plane 
exactly  at  loctd  dawn,  of  course,  due  to  the  tendency  of  outer-magnetospheric  field  lines  to 
be  swept  back  antisunward.  The  parameter  bound,  defined  in  a  similar  way. 

We  have  made  the  physical  assumption  that  the  equatorward  edge  of  the  electric-field- 
reversal  region  of  the  ionosphere  maps  to  the  equatorial  locations  specified  in  equation 
(2.4-2).  In  fact,  it  is  not  clear  observationally  how  far  out  in  the  tail  the  equatorward  edge 
of  the  electric-field-reversal  region  maps.  Although  our  choice  is  sensible,  it  may  not  be 
very  accurate. 

In  any  case,  we  assume  that  boundary  b  (ellipse  2)  is  an  ellipse  in  the  ionosphere  in 
flat-polar  coordinates,  and  we  determine  the  exact  location  of  that  ellipse  by  using  the 
magnetic-field  model  to  map  from  the  equatorial  plane  to  the  ionosphere.  Note  that  the 
ionospheric  location  of  ellipse  2  is  defined  by  magnetic-field  mapping  from  the  equatorial 
plane,  not  from  a  measured  boundary. 

Our  magnetic-field-model-based  algorithm  for  locating  ellipse  2  was  chosen  for  two 
practical  reasons: 

(i)  The  input  data  stream  that  is  presently  expected  to  be  available  for  operational  use  of 
the  MSM  does  not  include  the  location  of  the  equatorward  edge  of  the  electric-field-reversal 
region  or  (roughly  equivalently)  the  equatorward  edge  of  the  region- 1  Birkeland  currents. 
Routine  automatic  determination  of  these  boundaries  is  somewhat  difficult. 

(ii)  There  is  always  a  chance  that,  for  operations  in  extreme  magnetospheric  conditions  or 
with  sparse  or  bad  input  data,  our  magnetic-field  model  will  be  highly  inaccurate.  If  we 
chose  our  main  modeling  boundary  from  an  observed  ionospheric  boundary,  the  equatorial 
map  of  that  boundary  might  be  ridiculous  from  the  viewpoint  of  particle  tracing:  either  on 
field  lines  that  extend  to  unphysically  large  distances  from  Earth  or  field  lines  that  lie  inside 
synchronous  orbit.  Either  of  these  circumstances  would  prevent  the  particle  tracer  from 
functioning  and  would  compromise  the  primary  objective  of  the  model.  Our  conservative 
algorithm  for  locating  ellipse  2  essentially  eliminates  the  probability  of  this  type  of  failure. 

We  have  found  that  our  magnetic-field-model-based  algorithm  for  locating  boundary 
b  sometimes  leads  to  ionospheric-electric-field  patterns  that  are  not  optimally  reahstic.  We 
have,  therefore,  developed  an  algorithm  for  post-correcting  the  ionospheric-electric-field 
pattern,  after  the  magnetospheric  particle  traces  have  all  been  completed,  using  the  location 
of  the  main  ionospheric  electric-field  reversal  as  additional  input  information.  Since  this 
information  is  not  expected  to  be  included  in  the  data  stream  that  will  be  routinely  available 
to  the  MSM,  we  have  not  included  that  correction  algorithm  in  the  MSM  package. 
However,  it  is  available  on  request. 

2.42.2.  Location  of  Ellipse  3,  which  is  Boundary  c 

To  estimate  this  boundary  from  the  available  observations,  we  make  use  of  its  close 
proximity  to  the  equatorward  edge  of  the  diffuse  aurora,  a  boundary  that  has  been  studied 
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extensively.  From  Figure  6  of  the  paper  by  Gussenhoven  et  al.  [1983],  we  estimate  that 
the  dawn  equatorward  edge  is  about  1°  poleward  of  the  midnight  equatorward  edge,  on  the 
average.  The  duskside  equatorward  edge  lies  several  degrees  further  poleward.  However, 
we  expect  the  electric-field  boundary  to  lie  equatorward  of  the  boundary  of  the  electron 
precipitation,  because  ions,  which  do  most  of  the  shielding,  penetrate  better  on  the  dusk 
side.  Therefore,  for  the  purpose  of  setting  the  equatorward  boundary  of  the  sunward-flow 
region,  we  choose  the  equatorward  edge  at  dusk  also  to  be  1°  poleward  of  the  midnight 
equatorward  edge,  i.e.. 


Table  2  of  Gussenhoven  et.  al.  (1983)  gives  straight-line  Kp  fits  to  the  noon  and  midnight 
equatorward  edges,  and  we  use  those  to  estimate  the  noon-midnight  difference.  We  define 
"noon"  to  be  the  average  between  1 100-1200  and  1200-1300,  and  "midnight"  to  be  the 
average  between  2300-2400  and  0000-0100.  We  eliminate  Kp  from  the  fit  formulas  and 
write  the  noon  equatorward  edge  simply  as  a  function  of  the  midnight  equatorward  edge. 
The  resulting  formula  is 


.eq.edge  _  ,eq.edge  J66.95-/C^i%)  X  1 .725  +  .SO 

t^noon  —  '\nidnight  1015 


Simply  combining  (2.4-l)-(2.4-4)  does  not  offer  a  positive  guarantee  that  ellipse  2 
lies  poleward  of  ellipse  3,  since  ellipse  2  is  derived  from  the  magnetic-field  model  and 
ellipse  3  is  derived  from  real-time  data  and  statistical  rules.  On  the  other  hand,  the  electric- 
field  model  will  do  something  ridiculous  if  ellipse  3  slips  poleward  of  ellipse  2.  Ridiculous 
results  are  also  guaranteed  if  the  two  boundaries  get  unphysically  close  to  each  other. 
Therefore,  we  enforce  minimum  thicknesses  for  region  2,  for  noon,  dawn,  dusk,  and 
midnight.  First,  we  set  some  (arbitrary)  input  parameters,  namely 


A/C„  =  A^drugHi  =  3.0°  (2.4-5) 


A  a'^  a  a"^ 

^^usk  ~  ^^awn 


45x10'^  PCP 


(2.4-6a) 


where  PCP  is  the  polar-cap  potential  drop  in  volts,  and  Emax,  the  maximum  allowable 
value  of  the  average  dawnward  electric  field  across  region  2  at  dawn  or  dusk,  is  given 
(arbitrarily)  by 


Emax  =  0.1  V/m 


We  now  define 


5,  = 


.  mod. bound.  .  .min.  .eg. edge 

A,  -  A  A,  -  A, 


for  i  =  noon,  dawn,  dusk,  and  midnight.  In  each  case,  if  6,  <  0,  we  then  set 


j^eg.edge 


A. 


eq.edge 


+  6, 


(2.4-6/?) 


(2.4-7) 


(2.4-8) 
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2.42  J.  Calculation  of  A,  B,  DX,  DY  for  Ellipses  2  and  3. 

Having  thus  located  the  latitudes  at  which  ellipses  2  and  3  cross  the  noon,  dusk, 
midnight,  and  dawn  meridians,  we  have  the  information  needed  to  calculate  the  parameters 
A,  B,  DX,  and  DY  for  those  ellipses.  The  formulas  used  are  the  following: 


DXi^L)  —  -  ^  [■^noon  ■ 

(2.4-9a) 

DY{D  =-  ^[f^usk- 

(2.4-96) 

X^  =  ^[l80°  - 

(2. 4- 10a) 

Y^  =  ^[m°-/^awn-/iusk] 

(2.4-106) 

A(L) 


BiD 


j  (y  -  [DX{L)]^  [DYjL)]^ 

V  [Z)y(L)]2 

/  {Xh'^iY^)^-[DX{L)]^[DY{L)\^ 

V  [DX(L)]2 


(2.4-llfl) 


(2.4-IIZj) 


2. 4. 2. 4.  Calculation  of  A,  B,  DX,  DY  for  Ellipse  1. 

As  discussed  in  section  2.4.3,  our  treatment  of  the  potential  distribution  in  the  region 
poleward  of  our  main  modeling  region  is  based  on  a  computerized  version  Rich  and 
Maynard  [1989]  of  the  Heppner  and  Maynard  [1987]  empirical  model  of  the  ionospheric 
potential  distribution.  To  utilize  that  model  within  our  formalism,  we  had  to  draw  ellipses 
on  each  of  the  seven  Heppner-Maynard-Rich  patterns  that  correspond  to  our  ellipses  1  and 
2.  From  those  ellipses,  the  drawing  of  which  involved  some  subjectivity,  we  derived 
parameters  A,  B,  DX,  and  DY  ,  which  we  label,  in  the  program  AHM,  BHM,  DXHM,  and 
DYHM.  These  parameters,  which  are  given  in  Table  2.4- 1 ,  have  two  indices,  the  first 
indicating  ellipse  number,  and  the  second,  which  is  called  IP  ATT,  running  from  1  to  7  and 
indicating  the  Heppner-Maynard  pattern  number.  These  Heppner-Maynard-Rich  patterns 
are  scaled  to  fit  into  our  model  polar  cap,  as  discussed  in  Section  2.4.3.  We  estimate  the 
width  of  the  electric-field-revers^  by  assuming  that  that  width  scales  in  the  same  way  as  the 
entire  polar  cap  does.  The  resulting  formulas  for  calculating  the  parameters  for  ellipse  #1 
are  as  follows 


>1(1)  =  >1(2) 


^AHM(\,IPATT) 

AHM{2,IPATT) 


B{\)  =  5(2) 


^BHMjlJPATT) 

BHM{1,IPATT) 


(2.4- 12a) 


(2.4-126) 
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DX(l)  =  DX{1)+ - - [DXHM{\,IPATT)-DXHM{2,IPATD] 

AHMi2,IPATT) 

(2.4- 12c) 


DYi\)  =  DYi2)+ - ^ - [Dr//M(l,/PA7T)-Dy//M(2,/F/47T)] 

BHM(2,IPATT) 

(2.4-1 2d) 


Table  2.4-1.  Heppner-Maynard  Ellipse  Parameters 


Ellipse 

IPATT 

AHM 

BHM 

DXHM 

DYHM 

1 

1 

17.45 

14.26 

-2.66 

1.60 

1 

2 

12.13 

13.78 

-5.53 

0.80 

1 

3 

16.06 

14.31 

-3.19 

1.12 

1 

4 

13.72 

13.78 

-0.11 

0.05 

1 

5 

14.79 

14.73 

-2.45 

0.48 

1 

6 

12.82 

11.70 

-1.65 

3.40 

1 

7 

15.00 

12.07 

-2.23 

1.22 

2 

1 

20.00 

16.97 

-3.30 

1.33 

2 

2 

16.17 

16.70 

-7.45 

0.32 

2 

3 

18.88 

17.07 

-4.10 

1.54 

2 

4 

19.31 

17.23 

-1.76 

0.85 

2 

5 

18.30 

17.93 

-2.77 

0.69 

2 

6 

16.60 

17.50 

-2.98 

1.86 

2 

7 

17.93 

16.60 

-3.14 

0.21 

2.4.3.  Potential  in  Region  0  (Polar  Cap). 

We  directly  used  the  subroutine  EPOT  that  was  kindly  supplied  by  F.  J.  Rich  for  this 
region.  The  mathematical  techniques  used  in  that  subroutine  were  described  by  Rich  and 
Maynard  [1989],  and  will  not  be  discussed  here.  We  just  describe  here  our  procedure  for 
scaling  and  sliding  the  Heppner-Maynard-Rich  patterns. 


Suppose  that  the  H-M-R  model  implies  that  boundary  a  satisfies  the  equation 

(xHM-DXHMf  ^  {yHM-DYHMf  ^  ^ 

A Bhm 

and  our  real-time  obsen-ational  data  indicate  that  boundary  a  should  satisfy 

{x-DX)^  ^  {y-DYf  ^  j 
A^  B^ 


(2.4.13) 


(2.4.14) 


Here  xhm,  ynM,  y,  are  related  to  the  corresponding  colatitude  0  and  local  time  <))  by 
relations  like  (lb)  anci  (Ic).  We  define  the  "rule  of  corresponding  points"  to  be 
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XHM  -  DXHM  ^  x-DX 
A-hm  ^ 

yHM  -  DYHM  ^  y-DY 
Bhm  B 


(2.4.15a) 


(2.4.15ft) 


At  the  beginning  of  our  modified  version  of  the  Heppner-Maynard-Rich  subroutine 
(called  EPOT),  we  do  a  change  of  coordinates.  Given  the  x  and  y  for  which  we  would  like 
to  compute  a  potential,  we  use  ^uation  (2.4.15)  to  compute  corresponding  and  yf^/^ 
values  for  the  corresponding  points  in  the  Heppner-Maynard-Rich  pattern.  From  Xf^^  and 
YnM’  wc  compute  the  corrected  geomagnetic  latitude  and  local  time  of  the  corresponding 
point  in  the  H-M-R  pattern,  and  we  use  those  values  in  the  main  H-M-R  calculation.  We 
then  use  the  original  EPOT  algorithm  to  calculate  the  potential  Vhm  at  point  (jc//a/  .  Ynm)- 
A  final  scaling  process  scales  the  potentials  themselves  for  consistency  with  the  polar-cap 
potential  drop  PCP  that  has  been  estimated  on  the  basis  of  real-time  input  data: 


V 


Vhm 


_ PCP _ 

VMAXilPATT)  -  VMINUPATT) 


(2.4.16) 


where  VMAX  and  VMIN  are  the  maximum  and  minimum  potential  values  in  the  Heppner- 
Maynard-Rich  pattern  for  the  relevant  pattern.  The  values  of  VMAX  and  VMIN  ,  which 
we  read  off  matrices  representing  the  Heppner-Maynard  patterns,  are  given  in  Table  2.4-2. 


Table  2.4-2.  Maximum  and  Minimum  Potentials  in  the  Heppner-Maynard-Rich  Patterns 


IPATT 

VMAX 

VMIN 

1 

-42280 

2 

55354 

-16003 

3 

14390 

4 

11287 

TtrriM 

5 

9329 

-12947 

6 

13249 

-14428 

7 

™lll« 

-13460 

2.4.4.  Potential  in  Region  3  (Middle  and  Low  Latitudes). 

2. 4. 4.0.  Introductory  Comments 

A  full  and  realistic  model  of  the  subauroral  electric  field  would  include  two 
components:  the  quiet-time  field  and  the  disturbance  field.  The  standard  representation  of 
the  quiet-time  field  is  an  observational  model  developed  by  Richmond  et  al.  [1980],  based 
on  quiet-time-average  incoherent-backscatter-radar  data.  We  have  not  included  this  quiet¬ 
time  electric-field  distribution,  because,  being  a  static  field,  it  has  essentially  no  effect  on 
the  transport  of  magnetospheric  particles,  which  is  our  primary  objective. 
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2. 4. 4.1.  General  Expression  for  the  Low-Latitude  Potential. 

The  disturbance  field,  on  the  other  hand,  has  an  important  effect  on  the  transport  of 
magnetospheric  particles,  and  must  be  included  in  the  MSM.  Unfortunately,  the 
disturbance  field  is  not  really  very  well  determined  or  understood  at  this  point,  and  we 
expect  this  to  be  a  major  focus  of  research  in  the  next  few  years,  research  involving 
collaborations  between  thermosphere  and  ionosphere  modelers  and  our  group.  At  the 
present  time,  we  think  that  we  probably  have  a  reasonable  idea  of  the  electric-field  patterns 
corresponding  to  prompt  penetration  of  magnetospheric  effects  to  low  latitudes  [Spiro  et 
al.,  1988].  Those  patterns  can  be  roughly  parametrized  as  follows: 

oe 

Vp3(0.<l>)  =  <V’(r)>  +  F(r)  sin-P(8)  [flm  sin(m(l))  +  cos(m(}))]  (2.4-17) 

m=0 


where  F(t)  represents  the  overall  strength  of  the  electric-field  penetration  at  time  t. 


To  estimate  the  a^s  and  b^s,  we  utilized  a  least-squares-fit  procuedure  to  match  the 
equatorial  electric  fields  from  a  run  carried  out  with  the  Rice  Convection  Model,  specifically 
run  3  from  the  paper  of  Spiro  et  al.  [1988].  The  result  was 


Mt’U  ^  0-5249  cos  (0)  -  0.0265  cos (20)  -  0.0541  cos  (30)  + 
-1-  0.0939  sin(0)  -  02883  sin(20)  -i-  0.0810  sin(30) 


(2.4-18) 


where  the  eastward  electric  field  £,j,  is  in  mV/m,  and  0  is  our  usual  local-time  angle,  0  at 
noon,  etc.  The  maximum  westward  electric  field  implied  by  (2.4-18)  is  0.86  mV/m. 
Pulling  the  maximum  westward  electric  field  out  of  the  formula  gives 


£<)  =  £w«t.max[ 0.6 103  cos(0)  -  0.0308  cos(20)  -  0.0629  cos(30)  -i- 
+  0.1092  sin(0)  -  0.3352  sin(20)  -h  0.0942  sin(30) 

Figure  2.4-2  shows  /  E^est,max  a  function  of  0,  along  with  the  integral 


(2.4-19) 


1 


£(tid0 


'west  .max 


=  0.6103  sin(0)  -  0.0154  sin(20)  -  0.0210  sin(30)  -  0.1092  cos(0)  -i- 


+  0.1676  cos{2<^  -  0.0314  cos(3(P) 


(2.4-20) 


Note  that  the  normalized  electric  field  is  strongest  in  the  region  between  midnight  and 
dawn,  where  it  is  negative  (westward).  The  total  variation  of  the  normalized  integral,  the 
difference  between  its  maximum  and  minimum  values,  is  given  by 


E 


wesi.max 


128 


(2.4-21) 
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Expressions  (2.4-20)  and  (2.4-21)  are  assumed  to  hold  at  any  latitude  equatorward  of  the 
shielding  layer,  because  we  are  assuming  separability  of  the  0  and  (j)  dependences.  Our 
comput^  RCM  patterns  exhibited  a  rough  separability. 

The  potential  along  the  shielding  layer  is  then  given  approximately  by 

=  -sin(e,/u.w)/?/|  d(I)' £:0((t)’)  (2.4-22) 

where  ^s}udd  is  the  colatitude  of  the  shielding  layer  and  /?/  is  the  geocentric  distance  of  the 
ionosphere  (®65(X)  km). 


Substituting  (2.4-20)  in  (2.4-22),  and  assuming  the  same  power-law  dependence  on 
sin(0)  as  before,  we  obtain 


S  in(  Q  shield) 

S  in(  ^shield) 

sin(25°) 

sin(0) 

yiZ^i^shUidA)  =  -  2.75 


X  [0.6103  sin(<|))  -  0.0154  sin(2<{))  -  0.0210  sin(3(l>)  -  0.1092  cos((t))  +  (2.4-23) 
+  0.1676  cos(2(t)) -0.0314  cos(3(|))  ]  +  VBAR 


Setting  p  =  1.38  gave  the  best  fit  to  the  computed  penetration  computed  in  the  SUNDIAL 
RCM  runs.  In  ^e  program,  we  have  neglected  the  dependence  of  ^shield  on  and  let 
^shield  =  A(3)-DX(3),  which  is  actually  appropriate  for  local  midnight. 

2.4.42,,  EstiffUltiOK  of  Eyvest,  max{ ^shield) ■ 

We  use  the  motion  of  the  equatorward  boundary  of  auroral-zone  electron  precipitation 
as  an  indicator  of  the  direct  electrostatic  penetration  of  magnetospheric  electric  fields 
through  the  shielding  layer. 

Assume  that  the  equatorward  edge  of  the  auroral  electrons  is  a  convection  boundary  -- 
the  inner  boundary  of  the  plasma-sheet  electrons  that  have  recently  come  in  from  the  tail. 
Now  consider  the  motion  of  the  ionospheric  mapping  point  of  electrons  that  lie  on  that 
boundary.  The  motion  of  these  boundary-electrons  is  a  combination  of  gradient/curvature 
drift  (mostly  eastward)  and  ExB-drift.  The  latitudinal  motion  of  the  ionospheric  map  of 
these  boundary  electrons  is,  we  assume,  due  predominantly  to  £xfi-drift  that  results  from 
the  east- west  component  of  the  ionospheric  electric  field: 

RjdA  ^  .^wesL  (2.4-24) 

dt  Bir 

where  /?/  is  the  geocentiic  radius  of  the  ionosphere  (6500  km,  say),  A  =  latitude  in  radians 
of  the  northern-ionosphere  crossing  point  of  the  boundary,  and  is  the  absolute  value  of 
the  radial  component  of  the  magnetic  field,  also  at  ionospheric  altitude.  Setting  = 
.5x10'^  T,  (2.4-24)  can  be  rewritten  as 

E^eJf'  =  -1576  (2.4-25) 

dthr 


where  A"  is  now  latitude  in  degrees. 
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If  we  use  (2.4-25)  to  estimate  of  the  maximum  westward  electric  field  at  the  shielding 
layer  --  specifically,  the  equatorward  edge  of  the  shielding-layer  region,  we  are  being 
sloppy  in  several  respects,  notably  the  following: 

1.  The  electron  inner  edge  is  not  necessarily  at  the  equatorward  edge  of  the  shielding- 
layer  region.  On  the  dusk  side,  there  is  frequently  a  gap  between  the  equatorward  edge  of 
the  electrons  and  the  equatorward  edge  of  the  shielding;  this  is  region  of  the  SAID 
(Subauroral  Ion  Drift)  events  [Spiro  et  ai,  1978].  By  equating  the  electron  inner  edge 
with  the  shielding  layer,  we  are  effectively  neglecting  the  potential  drop  across  the  rapid- 
trough-flow  region  compared  to  the  potential  drop  along  the  shielding  layer.  This  error  will 
tend  to  make  us  overestimate  the  penetration  effect. 

2.  If  precipitation  becomes  so  strong  that  the  electron  equatorward  edge  is  essentially  a 
precipitation  boundary,  rather  than  a  convection  boundary,  then  (2.4-25)  will  probably  tend 
to  underestimate  the  penetration  electric  field. 

3.  Our  procedure  assumes  implicitly  that  the  meridional  motion  of  the  inner  edge  can  be 
equated  to  the  meridional  motion  of  an  individual  inner-edge  particle;  this  is  valid  if  the 
meridional  motion  occurs  on  a  time  scale  that  is  small  compart  to  the  particle  drift  time; 
otherwise,  the  particle's  equatorward  motion,  for  example,  may  be  limited  by  the  length  of 
time  that  it  is  subjected  to  the  nightside  westward  electric  field.  In  this  respect,  our  simple 
approximation  is  best  when  the  plasma-sheet  inner  edge  is  moving  fast,  and  penetration  is 
strong.  Of  course,  those  are  the  conditions  under  which  it  is  important  to  get  a  good 
estimate  of  the  penetration. 

Substituting  (2.4-26)  in  (2.4-23)  gives  an  explicit  expression  for  the  low-latitude 
potential.  Using  the  fact  that  the  square-bracketed  quantity  in  (2.4-23)  has  a  total  range  of 
1.28,  we  find  that  the  total  potential  drop  across  the  shielding  layer  (the  "peneu'ation 
potential")  is  given  approximately  by 

v!^^enet  =  *  13.1  g  ^SmOsHield)  (2.4-26) 

dthr 


where  Q^huid »  which  is  an  effective  shielding-layer  colatitude,  is  taken  to  be  A(3)-DX(3). 
Note  that  is  defined  to  be  positive  in  times  of  increasing  convection,  i.e.,  auroral 
zone  moving  equatorward. 

2.4.43.  Estimation  of  VBAR. 

It  is  well  established  that  the  there  is  more  sunward  flow  on  the  duskside  auroral 
zone,  on  the  average,  than  on  the  dawn  side.  This  flow  occurs  in  what  we  call  region  2, 
plus  the  sunward-flow  part  of  region  1.  Recent  work  by  Lu  et  al.  [1989]  suggests  that 
ratio 


Potenad  drop  across  duskside  sunward-flow  region 
Potential  drop  across  downside  sunward-flow  region 


is  about  1.5,  on  the  average.  We  choose  the  value  of  VBAR,  the  average  potential  at  low 
latitudes,  such  that  the  ratio  of  the  potential  drops,  in  the  auroral  zone,  is  held  at  a  value  of 
1.5.  The  potential  drop  across  the  auroral  zone  is  approximately  VBAR  -  0.38  Vpg„g,  - 
on  the  dusk  side  and  approximately  V,^  -  VBAR  -  0.62V on  the  dawn  side. 
Requiring  that  the  ratio  of  these  two  potential  drops  be  1.5  yields  the  following  expression 
for  VBAR: 


penet 


VBAR  =  0.6  V,^  +  0.4  V -  022  V 


(2.4-27) 
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2.4.5.  Potential  in  Region  2  (Main  Sunward  Flow  Region  of  the  Auroral 
Ionosphere). 

Because  the  low-latitude  field  is  normally  rather  weak,  but  it  is  crucial  to  keep 
appropriate  continuity  across  the  shielding  layer,  we  simply  add  an  extra  (usually  large) 
aurora-zone  field  to  a  smooth  extrapolation  of  the  potential  used  for  region  3: 


V2(0,<{»)  =  +  F«(e,(t))  (2.4-28a) 


where  V/ow.*  is  the  extrapolation  into  the  auroral  zone  of  the  low-latitude  potential,  and  Vaz 
is  the  main  potential  of  the  auroral  zone.  We  carry  out  the  smooth  extrapolation  in  the 
following  simple  way: 

V/„w,x(0.<{»)  =  +  (0-0c)  (2.4-28/7) 

30c 

For  the  main  auroral-zone  potential,  we  use  the  formula 


1-1-1- 


(0.-0) 


2\r 


V«(0.(|))  =  V<,,«((j))i 


+  Fc„,(e,(|))  (2.4-29a) 


1 


46' 


(J 


where,  in  the  first  term  on  the  right  side,Vi,  a2  is  the  auroral-zone  contribution  to  the 
potential  at  boundary  b,  and  we  choose  the  exjwnent  r  to  be  equal  to  unity.  In  the  second 
term,  which  is  intended  to  represent  the  rotation  of  equipotentials  in  the  Harang- 
discontinuity  region. 


=  .  6.75  4  W 


where 


1 7i(<l>-7t)  \ 

<|>  =  cos -  if 

^idthl 


2(7t-(l)) 


^  ^idfh 

0'  =  0  otherwise 


iQc-%y 


<  1 


and  our  present  choice  for  the  A0  coefficients  is  as  follows: 

A0M-irft/t  =  10 

=  2.0 

We  take  r  =  1  and 


A0  =  (0c  -  Qb) 


1.25  -  0.75  CO 


(2.4-29/7) 


(2.4-29C) 


(2.4-29d) 


(2.4-30) 


which  means  that  A0  =  (0c  -  0i,)/2  for  0  =  .757t  and  A0  =  2  (0c  -  0/,)  for  0  =  1.757c. 
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Figure  2.4-2  shows  some  typical  electric-field  curves  produced  by  this  functional 
form.  Specifically,  the  curves  shown  give  meridional  electric  fields,  with  the factor 
set  equal  to  unity.  Also,  <()'  is  set  to  zero. 

2.4.6.  Potential  in  Region  1  (The  Electric-Field-Reversal  Region). 

We  treat  region  1  as  a  kind  of  transition  region:  we  design  the  potential  function  for 
region  1  so  that  it  fits  smoothly  onto  both  region  0  and  region  2.  The  procedure  is  the 
following: 

(i)  Find  0a(<I))  and  0i,((())  for  the  <|)-value  in  question; 

(ii)  Find  Vq  =  VqI ej(t)),<l>],  Vq'  =  dV qI 6.0]/dOlea((l))>  ^2  =  ^2^ ^bW><Ph  and 
V2' s  dV2(d,^)/ddlgf^^)  using  the  analytic  formulas  described  for  regions  0  and  2; 

(iii)  Set 

V(0,<|>)  =  no+  oo'  (0-0«)  +[3(u2-do)-(2uo'+^^2’)80] 

3  S0  (2.4-31) 

-i-[2(i;o-t'2)+(t>o’+»2’)50] 

50^ 

where  50  s  0|,(<j))  -  0a(<t>) .  This  formula  makes  both  V  and  its  derivative  continuous  at 
both  0a(<|))  and  0i,(<|)). 

2.4.7.  Potentials  and  Electric  Fields  on  the  Boundaries. 

It  is  useful  to  start  the  potential-calculation  procedure  by  calculating  the  potentials  and 
their  normal  derivatives  on  boundaries  a,  b,  and  c,  and  storing  those  values,  before 
performing  the  main  potential  calculations. 

For  boundary  a ,  the  poleward  edge  of  the  electric-field-reversal  region,  we  calculate 
the  potential  directly  by  calling  EPOT,  the  subroutine  that  calculates  V  in  region  0,  for 
latitudes  and  longitudes  on  the  boundary.  To  calculate  the  normal  derivative  we  call  the 
same  subroutine  a  second  time,  this  time  for  points  a  latitudinal  distance  e  inside  the 
boundary.  The  normal  derivative  is  computed  by  subtracting  the  two  potentials  £  apart, 
and  divi^ng  by  e. 

For  boundary  b,  equatorwaid  edge  of  the  electric-field-reversal  region,  we  calculate  V 
by  calling  EPOT  again,  this  time  for  colatitude  0^,. 

For  boundary  c,  the  equatorwaid  edge  of  the  shielding  layer,  we  calculate  both  V  and 
its  normal  derivative  by  calling  the  low-latitude  routine,  for  points  on  the  boundary  and  also 
for  points  a  distance  £  equatorwaid  of  boundary  c. 

This  provides  sufficient  information  to  completely  define  the  potential  in  region  2. 
That  procedure  guarantees  continuity  of  the  potential  and  its  derivative  at  boundary  c,  and 
of  the  potential  at  boundary  b.  By  remembering  the  previously  calculated  potential  on  the 
boundary  and  calling  the  region-2  routine  for  a  point  £  equatorward  of  boundary  b,  we 
establish  the  normal  derivative  of  V  at  that  boundary. 

Since  the  potential  and  its  normal  derivative  are  now  defined  at  both  boundaries  a  and 
b,  there  is  sufficient  information  available  for  complete  determination  of  the  coefficients  in 
the  latitude  expansion  of  the  expansion  in  region  1  (equation  (2.4-31)). 
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2.4.8.  Numerical  Illustrations. 

Figures  2.4-3  to  2.4-5  show  equipotential  diagrams  for  several  cases,  one  in  which 
the  equatorward  edge  of  the  auroral  zone  is  not  moving  and  there  is  thus  no  penetration  to 
low  latitudes  (Figure  2.4-3),  one  in  which  the  equatorward  edge  is  moving  equatorward 
(Figure  2.4-4)  and  there  is  consequently  a  westward  electric  field  across  the  night  side  at 
low  latitudes,  and  one  in  which  the  equatorward  edge  is  moving  poleward  (Figure  2.4-5) 
and  there  is  an  eastward  electric  field  at  low  latitudes. 
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2.5.  Initial-Condition,  Boundary-Condition,  and  Reference  Fluxes 
2.5.0.  Introductory  Comments. 

For  the  purpose  of  setting  boundary  and  initial  conditions  for  the  MSM,  and  for 
setting  upper  and  lower  limits  on  the  fluxes,  we  have  calculated  a  flux  matrix  called 
FLXMAT,  which  is  based  loosely  on  statistical  studies  and  previously  published 
observational  data.  We  first  discuss  the  construction  of  the  FLXMAT  matrix,  then  discuss 
the  way  in  which  we  use  it  for  initial  conditions,  etc. 

There  is  a  FLXMAT  matrix  for  each  of  three  chemical  species  (electrons,  H+,  O"^). 
Each  of  these  matrices  specifies  values  of  the  differential  flux  (cm-2  sr^  eV-i)  for  seven 
Kp  values  (0  to  6),  for  29  energy  values  from  10^  eV  to  10^  eV  in  steps  of  100-25,  and  for 
the  following  four  geocentric  distances: 

(i)  r  =  13,  which  is  supposed  to  represent  a  full- strength  plasma  sheet. 

(ii)  r  =  6.6,  which  represents  the  synchronous-orbit  region. 

(iii)  r  =  4,  which  represents  (roughly)  the  peak  of  the  outer  radiation  belt. 

(iv)  r  =  3,  which  represents  (roughly)  the  slot  region  between  the  inner  and  outer  belts. 

In  all  cases,  the  energy  dependence  of  FLXMAT  is  assumed  to  take  the  form  of  a 
kappa  distribution.  We  begin  the  discussion  by  reviewing  the  properties  of  the  kappa 
function.  The  coefficient  values  assumed  for  FLXMAT  for  each  of  the  four  values  of  r 
will  then  be  presented  in  later  sub-sections. 

2.5.1.  Properties  of  the  Kappa  Distribution. 

The  basic  properties  of  the  kappa  distribution  are  as  follows  (from  Vasyliunas 
[1968]): 


/(„)  =  JL  nic^i) 


wj  (TtK)^-^  r(K-^) 


1  + 


lK+1 


(2.5-1) 


where  N  =  number  density.  The  differential  flux  (particles/area/time/energy/solid  angle)  is 
given  by 

j(E)  =  N  Wq  r(K+l) 


2  Eo  (7tK)^/2  r(K-^) 


1+- 


x+1 


(2.5-2a) 


where,  of  course,  Eo  =  mwo^/2  and  E  =  mv^ll.  In  program  units,  the  same  equation 
becomes 


/(£)  =  (1.68  X  10*  keV  *  cm-2  s  '  ster')  F(0 


[1  keVl 

1/2 

[  N  1 

jL' 

Eo  \ 

-1  cm-3. 

^Eo 

K+l 

where 

r(K-n) 

k3/2  r(K-i) 


(2.5- 


F{k)  = 


(2.5-2C) 
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The  function  F{k)  is  a  very  slowly  varying  function  of  k;  some  values  are  given  in  Table 
2.5-1  below. 


Table  2.5-1.  The  function  fTK") 


K 

f(k) 

2 

0.7978846 

3 

0.86862694 

4 

0.9027036 

5 

0.92274588 

6 

0.9359421 

7 

0.94528801 

8 

0.95225411 

9 

0.95764671 

The  energy  density  is  given  by 


and,  correspondingly, 


u  =  NEo}^ 
2  K  -  i 
2 


k-2- 


Eo  =  I 


where  is  the  average  energy  of  a  single  particle. 
2.5.2,  Calculation  of  F  LX  MAT  for  Electrons. 


(2.5-3fl) 

(2.5-3^) 


2.52.1.  FLXMAT  for  Electrons  at  r  =  13. 

Our  number  density  and  temperature  estimates  are  based  on  Figure  5  of  the  paper  by 
Huang  and  Frank  [1986].  We  choose  values  consistent  with  z  =  0,  i.e.,  as  close  as 
possible  to  the  center  of  the  central  plasma  sheet.  We  associate  AE  <  200  with  Kp  =  1  and 
500  <  AE  <  1000  with  Kp  =  5.  On  this  basis,  we  choose  the  number  density  and  ion 
temperature  as  follows: 


plasma  sheet  _  0-^(^P~l)  i5-Kp) 

^  _  {\2keW)iKp-\)  +  (6.5keW){5-Kp) 

i  I  ^  - - - 


(2.5-4) 


(2.5-5a) 


According  to  Baumjohann  et  al.  [1989],  the  ion  and  electron  temperatures  are  highly 
correlated,  with 


^  =  7.8  (2.5-5^) 


Combining  (2.5-5a)  and  (2.5-5fr)  and  multiplying  by  1.5  to  convert  to  average  electron 
energy,  we  obtain 
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^avg.  plasma  sheet  _  (2.31  keV)(Ap-l)  +  (1.25  keV)  (5-Kp)  (2  5-5c) 

For  the  plasma  sheet,  we  use  a  kappa  distribution,  with  k  =  6.  That  distribution  uses 
the  energy  parameter  £<,.  which,  according  to  (2.5-3b),  is  half  the  average  energy.  We 
therefore  obtain 


Eo  =  (0.2885  keV)  (Kp-l)  +  (.15625  keV)  (S-Kp)  (2.5-5d) 

The  flux  is  given,  in  our  program  units,  by  substituting  (2.5-4)  and  (2.5-5d)  in  (2.5-2^?).  I 
suggest  that  we  make  a  genei^  conservative  practice  of  using  Kp  =  6  values  for  Kp  >  6. 

25.2.2.  FLXMAT  for  Electrons  at  r=6.6. 

We  use  the  sum  of  two  k  =  6  distributions  to  represent  the  situation  at  synchronous 
orbit,  one  for  Garrett's  low-energy  electrons  and  a  second  for  his  high-energy  electrons. 
The  values  of  £<,  and  number  density  N  are  given  in  the  following  table.  These  are  local¬ 
time  averages  derived  from  Garrett’s  latest  model  [Garrett  and  Gaudet,  1989].  The 
standard  reference  on  this  line  of  research  is  Garrett  er  a/.  [1981]. 

Table  2.5-2.  Garrett-Model  Electron  Fluxes  at  Synchronous  Orbit 


Ep 

MSIBSiSSSM 

mmssm 

wunmm 

mM&sm 

g.2(cV) 

0 

0.5697 

41.382 

■illlH  ^ 

0.1826* 

6502.9 

4877.2* 

1 

0.6034 

54.404 

40.803 

0.2261 

6443.2 

4832.4 

2 

0.6292 

64.893 

48.670 

0.2583 

6418.6 

4814.0 

3 

0.7002 

95.573 

71.680 

0.3418 

6397.1 

4797.8 

4 

0.8135 

148.159 

111.119 

0.4602 

6421.3 

4816.0 

5 

1.0347 

256.939 

192.704 

0.6446 

6525.7 

6 

1.4468 

456.949 

342.712 

0.8505 

6767  0 

5075.2 

*  These  should  be  overriden  by  other  values.  (See  following  text.) 


However,  based  on  the  tendency  of  the  Garrett-model  to  overestimate  the  >30  keV 
electron  fluxes  during  quiet  times  in  the  April  1988  event,  we  need  to  override  Garrett's 
ATp  =  0  values  for  Ne2,  Te2.  and  Eoi-  For  this  purpose,  we  found,  in  the  literature,  seven 
cases  where  D.  N.  Baker  and  collaborators  had  published  spectral  information  on  30-300 
keV  electrons  for  a  well-defined  substorm  period,  including  a  quiet  period  before  the 
substorm.  The  dates  and  references  are  given  in  Table  2.5-3. 


Table  2.5-3.  Sources  of  Quiet-Time  Data  on  >  40  keV  Geosynchronous  Electrons 


Date 

ESIESIH 

Reference 

3/22/79 

14 

[McPherron  and  Manka,  1985] 

9/8/77 

1 

[Baker  et  ai,  1978] 

4 

[Baker  et  aL  1984  ] 

7 

[Baker  et  aL  1984] 
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2/1/83 

[Baker  et  al,  1984] 

2/1/83 

[Baker  et  al,  1984] 

7/29/77 

3 

[Baker  et  al,  1982] 

(The  notations  "sp.019"  and  "sp.025"  refer  to  two  different  geosynchronous  spacecraft,  at 
different  local  times,  observing  the  same  event.  The  results  from  the  two  spacecraft 
measuring  the  same  event  from  different  locations  differed  by  about  as  much  as  the 
dispersion  among  the  different  events,  so  we  counted  them  as  independent  measurements.) 

Since  we  seek  an  observation-based  representation  of  quiet-time  fluxes  at 
synchronous  orbit,  we  checked  the  Dst  index  for  each  of  these  times  to  make  sure  that 
none  of  them  occurred  during  a  magnetic  storm.  Only  the  last  listed  event,  the  CDAW-2 
event  of  7/29/77,  occurred  during  a  magnetic  storm.  We  therefore  eliminated  the  7/29/77 
event  from  consideration. 

In  each  case,  we  scaled  the  flux  values  from  the  beginning  of  the  substorm  interval 
plotted,  representing  the  quiet  time  before  the  arrival  of  the  substorm  fluxes.  In  some 
cases,  there  was  a  dramatic  dropout  of  fluxes  within  an  hour  of  the  substorm  expansion 
phase:  our  measured  values  were  always  taken  well  before  that  dropout.  We  then 
averaged  the  six  measured  quiet-time  log(flux)  values  and  fitted  the  average  to  a  k=6 
distribution.  The  comparison  of  the  averaged  data  points  and  the  kappa  function  are 
displayed  in  Figure  2.5-1. 

The  resulting  kappa-function  parameters,  which  we  use  to  override  the  tabulated  hot- 
plasma  Garrett  values,  are  as  follows: 

=  4.926  x  lO'^  (2.5-6a) 

£o2(eV)  =  21,0(X)  (2.5-66) 

The  rationale  for  overriding  Garrett's  values  with  these  values,  which  are  derived  from  >30 
keV  electron  data  for  pre-substorm  conditions,  is  as  follows.  We  are  going  to  use  Kp=Q 
values  as  a  base  (equilibrium)  level  for  the  trapped  radiation.  We  suspect  that  Garrett's 
Kp^  values  may  have  been  contaminated  somewhat  by  substorm  effects.  I  should  also 
comment  that  Garrett's  older  classic  paper  on  the  statistics  of  geosynchronous  electrons 
[Garrett  et  al.,  1981]  indicated  a  much  lower  value  of  Ngi  for  Kp  =  0  than  is  indicated  in 
the  first  line  of  the  preceding  Table.  However,  in  overriding  the  straight  Garrett  values, 
we  may  be  sacrificing  accuracy  in  the  representation  of  electrons  in  the  1-30  keV  range  to 
improve  accuracy  in  the  >  30  keV  range. 

25.23.  FLXMAT  for  Electrons  at  r  =  3  and  4. 

Fluxes  for  r  =  3  and  4  were  scaled  from  Figures  5-45  and  5-46  of  Spjeldvik  and 
Rothwell  [1985],  representing  the  AE7-HI  and  AE7-LO  NASA  Radiation  Belt  Models. 
We  considered  four  energy  channels  (40-100  keV,  100-250  keV,  250-500  keV,  and  500- 
750  keV)  and  averaged  the  values  from  the  two  plots,  and  then  fit  the  result  with  a  kappa 
function.  Figure  2.5-2  compares  the  kappa-function  fit  with  the  scaled  data  for  both  L  =  3 
and  4. 

The  analytic  form  of  the  kappa-function  fits  are  as  follows: 
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j^\E)  =  (4.795  X  10^  cm-2  s'l  sr^  keV-^) 


h  EkeV 


+  ^keV 


65.8 


3.8 


10[(A;’-2)/6]  {2.5-la) 


j’^iE)  =  (1.106  X  10^  cm-2  s-i  s^i  keV'i)  — - \Qmp-2)i(>]  {2.5-lb) 


32  J 


For  L  =  3,  this  corresponds  to  k  =  2.8  and  Eo  =  23.5  keV.  For  L  =  4,  it  corresponds  to  k 
=  2and£^=  16keV. 

A  remark  is  needed  concerning  the  assumed  /(^-dependence.  The  NASA  Radiation- 
Belt  Models  do  not  contain  information  on  the  dependence  of  the  fluxes  on  the  level  of 
magnetic-activity,  and,  as  a  result,  we  have  to  do  something  arbitrary.  Scanning  plots  in 
two  old  references  {Owens  and  Frank ,  1968;  VanAllen,  1968],  we  found  that  the  range  of 
flux  variations  caused  by  hour-to-hour  variations  in  magnetic  activity  is  about  a  factor  of 
10.  Since  the  average  Kp  is  about  2,  the  average  values  that  were  read  off  the  graphs  from 
the  NASA  models  were  taken  to  correspond  to  Kp=2. 


2.5.3.  Calculation  of  F  LX  MAT  for  and  Ions. 

2.5.3.O.  Introductory  Comments. 

The  MSM  keeps  track  of  two  ion  species,  and  0+,  separately,  because  the  two 
species  have  substantially  different  loss  rates.  We  construct  two  ion  FLXMAT  matrices, 
one  for  and  one  for  O"^.  As  in  the  case  of  the  electron  FLXMAT  matrices,  each  matrix 
has  elements  for  seven  different  /Kp-values,  29  different  energies,  and  four  different  radial 
distances.  However,  in  the  case  of  ions,  we  must  also  seek  composition  information. 

In  some  cases,  we  will  assume  that  the  two  ion  species  have  the  same  £<,  but  different 
densities.  We  define  the  ratios 


We  thus  have 

no+ 


(2.5-8) 

«//+ 

N 

(2.5-9) 

l-t-rn 

Nr„ 

\+r„ 

(2.5-10) 

253.1.  Calculation  of  FLXMAT  for  r  =  13  for  Ions. 

Our  number  density  and  temperature  estimates  are  based  on  Figure  5  of  the  paper  by 
Huang  and  Frank  [1986].  We  choose  values  consistent  with  z  =  0,  i.e.,  as  close  as 
possible  to  the  center  of  the  central  plasma  sheet.  We  associate  AE  <  2(X)  with  Kp  =  1  and 
500  <  AE  <  KXX)  with  Kp  =  5.  On  this  basis,  we  choose  the  number  density  and  ion 
temperature  according  to  equations  (2.5,4)  and  (2.5.5a).  Multiplying  by  1.5  to  convert  to 
average  electron  energy,  we  obtain 
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^as>g.  plasma  sheet  ^  (18  keVX/Cp-l)  +  (9.75  keV)(5-/^p) 

4 


(2^-11) 


For  the  plasma  sheet,  we  use  a  kappa  distribution,  with  k  =  6.  That  distribution  uses 
the  energy  parameter  £<,,  which,  according  to  (2.5-3h),  is  half  the  average  energy.  We 
therefore  obtain 


Eo  =  (225  keV)  (Afp-l)  +  (121875  keV)  {5-Kp) 


From  Figure  5  of  Lennartsson  and  Sharp  [1982],  we  estimate 


r„ 


0.6  (/i:p-l)  + 0.06  (5-^:p) 


r„  =  0.06  for  Kp<  1 


(2.5-12) 


(2.5.13) 


(The  low-^Tp  cutoff  is  to  prevent  from  going  negative.)  Figure  6  of  Lennartsson  and 
Sharp  [  1984]  implies  that  the  average  energies  for  and  H+  in  the  plasma  sheet  are  nearly 
equal,  so  we  can  associate  Eo^*  and  Eo°*  with  the  Eg  value  derived  in  equation  (2.5-12). 

Figure  2.5-3  shows  plasma-sheet  and  0+  fluxes  computed  according  the 
prescription  given  in  this  section. 

2.532.  Calculation  of  FLXM AT  for  r  =  6.6  for  Ions. 

We  use  the  sum  of  two  k  =  6  distributions  to  represent  the  situation  at  synchronous 
orbit,  one  for  Garrett’s  low-energy  electrons  and  a  second  for  his  high-energy  electrons. 
The  values  of  E^  and  number  density  N  are  given  in  the  following  table.  These  are  local¬ 
time  averages  derived  from  the  model  of  Garrett  and  Gaudet  [1989]. 


Table  2.5-3.  Garrett-Model  Ion-Flux  Parameters  at  Synchronous  Orbit 


Kp 

MStBSSSSM 

Ti,(eV) 

iBnssssnii 

0 

0.491 

120.1 

180.2 

0.414 

16938.5 

25408 

1 

0.519 

172.3 

258.4 

0.411 

17891.6 

26837 

2 

0.540 

223.6 

335.4 

0.408 

18601.3 

27902 

3 

0.597 

403.2 

604.8 

0.400 

20482.0 

30723 

4 

0.687 

762.8 

1144.2 

0.385 

23314.8 

mssam 

5 

0.853 

1543.6 

2315.4 

0.349 

28539.4 

42809 

6 

1.111 

2788.3 

4182.4 

0.290 

38108.6 

57163 

For  the  low-energy  ions  (component  1),  we  base  our  composition  estimates  on 
Figure  5  of  Lennartsson  and  Sharp  [1982],  obtaining 


rni 


0.7  (KpA)  +  0.4  (5-Kp) 
4 


(25-14) 


Since  Garrett  models  fit  to  a  Maxwellian,  we  take  a  large  value  of  Kj,  namely 
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Ki  =  10 


(2^-15) 


For  the  high-energy  ions  (Garrett's  component  2),  we  use  the  paper  by  Gloeckler 
and  Hamilton  [1987].  Table  VI  of  that  paper  gives  values  of  /„  for  "Quiet"  and 
"Disturbed"  times,  which  we  associate  with  Kp  =1  and  5,  respectively.  This  assumption, 
combined  with  the  assumption  that  is  linear  in  Kp^  leads  to  the  formula 


=  0-154  0.767  (Ap-1)  ^ 

r„  =  0.1  M  for  Kp<  1 


(2.5-16) 


Figure  9  of  Gloeckler  and  Hamilton  [1987]  suggests  some  noticeable  differences 
between  the  energy  spectra  of  the  H+  and  0+  ions  at  L  =  6.6.  However,  the  differences  are 
fairly  subtle,  and ,  for  the  sake  of  simplicity,  we  use  the  same  £^2  and  K2  for  both  species, 
specifically  the  £^2  values  listed  in  Table  2.5-3  and 

K2  =  4  (25-17) 

The  Garrett-model  analysis  was  based  on  the  assumption  that  all  of  the  observed  ions 
were  H'^.  His  instrument  basically  measures  energy  and  flux.  The  energy  measurement  is 
independent  of  mass,  but  interpreting  the  flux  in  terms  of  number  density  and  energy 
involves  an  assumption  about  the  mass.  Interpreting  0+  as  being  H"*^  causes  a  factor-ofi 
four  underestimate  in  the  number  density: 

Ncarretl  ^  riH^  +  ^  (25-18) 

4 

The  actual  total  number  density  is  then  given  by 


H total  ~  HQgrrett 


(2.5-19) 


Table  2.5-4.  Corrected  Garrett  Number  Densities 
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2.533.  Calculation  of  FLXM AT  for  r  =  4  for  Ions. 

For  the  radiation  belts,  we  know  of  no  ^p-based  statistical  study  analogous  to  those 
of  Garrett  and  Gaudet  [1989]  or  Huang  and  Frank  [1986].  Our  FLXMAT  consequently 
has  to  be  based  on  some  individual  cases.  In  each  case,  we  performed  a  rough  digitization 
of  an  observed  j{E)  curve,  measuring  the  fluxes  at  intervals  of  100-25.  We  then  fitted  a  bi¬ 
kappa  distribution  to  the  resulting  curve.  Table  4  gives  the  resulting  bi-kappa  parameters. 
Figures  2.5-5  to  2.5-7  compare  our  fitted  curves  with  the  digitized  observational  ones,  for 
the  three  Kp  levels  considered  (0, 5+,  and  9-).  The  source  for  the  Kp=Q  curves  was  Figure 
9  of  Gloeckler  and  Hamilton  [1987]:  that  figure  represented  an  average  over  12  quiet 
periods,  for  L  =  3-5.  The  source  for  the  Kp  =  5+  curve  was  Figure  12a  of  Gloeckler  and 
Hamilton  [1987],  which  pertained  to  L=3.7-4.7  for  a  pass  during  a  magnetic  storm  on 
September  5, 1984.  The  Kp  value  for  the  time  of  the  pass  was  5+.  The  source  for  the  Kp 
=  9-  period  was  Figure  8fl  of  Hamilton  et  al.  [1988],  which  pertained  to  a  pass  during  the 
very  large  magnetic  storm  of  February  1986.  The  pass  in  question  occurr^  just  after  the 
minimum  in  Dst,  and  Kp  for  the  interval  was  9-. 


Table  2.5-5.  Bi-kappa  Parameters  for  r  =  4. 


Kp 

Ki 

A^i(cm-3) 

£«/(keV) 

K2 

£o2(keV) 

0 

H+ 

6 

0.1 

3 

10 

0.375 

200 

0 

0+ 

4 

0.43 

5 

- 

0 

- 

5+ 

H+ 

10 

0.2 

15 

6 

1.9 

40 

5+ 

0+ 

10 

OJfT 

5 

6 

1.5 

9- 

- 

0 

- 

2.3 

60 

Cfr 

' 

0 

- 

10 

5.0 

30 

Several  comments  are  needed: 

1 .  These  curves  are  all  based  on  data  from  the  CHEM  instrument  on  AMPTE,  which 
apparently  had  a  lower-energy  threshold  of  about  30  keV.  fluxes  were  usually  not 
shown  even  down  to  30  keV.  One  should  not  place  much  confidence  on  FLXMAT  below 
30  keV,  since  those  values  are  based  essentially  on  extrapolation  using  the  kappa  functions. 

2.  Gloeckler  and  Hamilton  [1987]  did  not  show  0+ fluxes  below  about  50  keV, 
although  the  H+  fluxes  were  shown  down  to  about  5  keV.  The  "data"  shown  for  0+ 
Figure  3  at  for  low  energies  are  simply  the  H+  fluxes:  I  just  guessed  that  the  0+  and  H+ 
fluxes  would  be  about  the  same. 

3.  In  general,  there  was  considerable  non-uniqueness  in  these  fits.  Only  in  the  case  of 
the  quiet-time  H+  flux  was  the  fit  improved  in  a  major  way  by  the  use  of  two  kappa 
functions  rather  than  one.  In  a  few  cases,  the  fit  didn't  seem  to  be  improved  at  all  by  the 
second  kappa  function,  and  I  set  the  density  of  one  or  the  other  equal  to  zero.  Because  of 
all  this  arbitrariness  in  the  use  of  the  second  kappa  function,  it  would  be  unwise  to  do  a  Kp 
-interpolation  on  individual  parameters  like  Ni,  Eo2,  etc.  It  would  be  better  to  interpolate 
the  log(flux)  values  themselves. 

4.  Gloeckler  and  Hamilton  [1987]  and  Hamilton  et  al.  [1988]  displayed  total  energy 
density  vs.  r  for  quiet  times  and  for  the  specific  events  that  we  are  using  as  benchmarks. 
These  observed  total  particle  energy  densities  covered  the  energy  range  from  30-315  keV. 
We  integrated  our  kappa  functions  over  the  same  range,  and  found  discrepancies  of  less 
than  0.13  in  the  log(base  10),  except  for  the  quiet- time  fluxes,  where  the  accuracy  was 
lower.  The  discrepancy  in  the  log  was  about  0.23  for  quiet-time  0+  and  0.45  in  the  case  of 
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quiet-time  O"^,  which  was  the  case  where  a  large  fraction  of  their  measured  data  points 
seemed  to  represent  only  a  small  fraction  of  the  total  flux  in  the  range  from  30  to  3 15  keV. 

I  suggest  that  FLXMAT  for  r  =  4  be  computed  by  summing  the  two  kappa  functions, 
for  Kp  =  0,  5+,  and  9-.  Values  of  FLXMAT  for  other  values  of  Kp  should  be  computed 
by  linear  interpolation  of  the  logarithms. 

25.3.4.  Calculation  of  FLXMAT  for  r  =  3  for  Ions. 

We  could  find  no  published  spectra  from  AMPTE  for  r  ~  3,  so  that  a  detailed  fitting 
exercise  like  that  carried  out  above  for  r  ~  4  is  not  possible.  However,  the  AMPTE 
investigators  have  published  curves  specifying  the  r-dependence  of  the  total  ion  energy 
density,  curves  that  extend  to  r  <  3.  Therefore,  we  use  the  same  table  of  bi-kappa 
parameters  as  was  employed  for  r  =  4,  but  with  all  of  the  densities  scaled  to  change  the 
energy  and  number  density  by  the  factors  suggested  by  the  observational  plots.  For  ^e  Kp 
=  0  and  Kp  =  5-i-  cases,  we  checked  the  corresponding  ratios  for  the  number  densities;  the 
resulting  ratios  were  nearly  the  same  as  the  ones  derived  from  the  energy-density  curves, 
which  implies  that  the  changes  in  spectral  form  must  not  be  major.  Table  2.5-6  gives  those 
factors. 


Table  2.5-6.  Ratios  of  Ion  Number  Densities  at  r  =  3  and  4 


Kp 

Species 

V(r=3)/Mr=4) 

Source 

0 

H+ 

0.32 

Fie.  10,  Gloeckler  and  Hamilton  [19871 

0 

0^ 

0.15 

Fie.  10,  Gloeckler  and  Hamilton  [19871 

5+ 

H+ 

0.32 

Fig.  11,  Gloeckler  and  Hamilton  [  19871 

5+ 

0+ 

0.10 

Fie.  11,  Gloeckler  and  Hamilton  [19871 

9- 

H+ 

2.7 

Fig.  7d,  Hamilton  et  al.  [19881 

9- 

0+ 

2.4 

Fig.  7d,  Hamilton  et  al.  [19881 

Normally,  the  ion  fluxes  drop  off  sharply  with  decreasing  r,  inside  r  =  4,  as  is 
evident  in  the  fu'st  four  rows  of  the  table.  However,  in  the  case  of  the  great  storm  of 
February  1986,  the  fresh  injection  of  particles  extended  inside  L=  3.  In  most  of  the  main 
phase,  the  peak  energy  density  occurr^  between  L=  2  and  3.  Of  course,  we  are  committed 
to  having  FLXMAT  based  on  L=  3,  4,  6.6,  and  13,  which  means  that  fluxes  for  L<3  are 
computed  by  linear  extrapolation  using  the  values  at  3  and  4.  Thus  for  very  large  Kp  our 
model  would  imply  that  ion  fluxes  increase  all  the  way  in  to  the  Earth,  were  it  not  for  our 
policy  of  using  the  Kp  =  6  values  for  all  Kp  higher  than  6.  A  consequence  of  that 
convention  will  be  that  our  FLXMAT  will  never  represent  the  fact  that,  in  a  very  large 
storm,  fresh  ring-current  ions  typically  penetrate  to  L  <  3. 

The  following  table  gives  the  bi-kappa  parameters  for  r  =  3. 
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Table  2.5-7.  Bi-kappa  Parameters  for  r  =  3. 


Kp 

Species 

Kl 

Eoj(keV) 

K2 

£«2(keV) 

H+ 

6 

1 

0.12 

200 

0+ 

4 

0.06 

5 

- 

0 

- 

5+ 

H+ 

10 

15 

6 

0.6 

40 

5+ 

O 

10 

5 

6 

0.15 

20 

H+ 

- 

0 

- 

10 

6.3 

60 

9- 

0+ 

- 

0 

- 

10 

11.8 

30 

2.5.4.  Initial-Condition  Fluxes 

25.4.1.  Calculation  of  the  value  of  FLXMAT  at  grid  points. 

A  primary  use  of  the  FLXMAT  (Kp,  r,  E)  matrix  is  in  establishing  the  initial- 
condition  flux  at  each  grid  point.  Define  the  symbol  jKpU,  J,  IE)  to  represent  the  value  of 
FLXMAT  that  is  assigned  to  a  given  grid  point  {IJ)  ana  a  given  invariant-energy  level  IE. 
We  calculate  jKpiF  J>  IE)  using  a  three-dimensional  interpolation,  the  three  dimensions 
being  Kp,  r,  and  E.  In  each  case,  a  basically  linear  interpolation  is  used,  but  there  are  a  few 
special  features  of  the  interpolation  routine  that  require  brief  discussion: 

(i)  The  interpolation  on  Kp  is  simple  and  linear  in  the  log  of  the  flux.  Kp  =  2+  is 
indicated,  in  the  code,  by  2.3333333,  2-  by  1.6666667.  Also,  we  use  the  Kp  =  6  flux- 
values  for  conditions  in  which  the  real  Kp  exceeds  6. 

(ii)  The  interpolation  on  r  is  carried  out  in  log-linear  fashion.  For  r  <  4,  we  use 
logidjKp(UJE))  = 

logioj  j|  X  logi(iFLXMAT{Kp,4,E))  +  logioj^j  x  \ogidFLXMAT(Kp,3,E)) 

logioji] 

(2.5-20aj 


For  4  <  r  <  6.6,  we  use 
\ogidjKp{U,IE))  = 

logiojyj  X  logio(FLXMA7(/(rp,6.6,£))  +  logioj^j  x  \ogx^FLXMAT{Kp,A,E)) 

logioj^) 

{1.5-lQb) 

For  6.6^  r<  13,  we  use 
logidi/KpdJJE))  = 
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logic 

('’"1 

\6.6j 

X  logio(FLXAfAr(/i:p,13,£))  +  logi, 

"(-1 

\rul 

^xlogi(iFLXMAT{Kp,6.6,E)) 

(2.5-20C) 


And  for  r  >  13,  we  use 

logio(/Arp(W£))  =  =  \ogiiFLKMAT{Kp,\XE))  +  \og^^]^  (2.5-20d) 

The  latter  condition  just  ensures  a  reasonable  decline  of  pressure  with  distance  out  into  the 
plasma  sheet.  Specifically,  it  has  the  density  and  pressure  declining  as  1/r,  with  the 
temperature  remaining  constant.  Considering  Ae  present  confusion  about  the  variation  of 
temperature  with  density  in  the  plasma  sheet,  this  seems  like  a  reasonable  simple  choice. 

(iiil  Interpolation  in  energy.  The  kinetic  energy  of  a  particle  of  species  IE  at  grid  point 
(IJf)  is  given  by 


£(/,  J,  IE)  =  \ALAM(IE)\  x  VM(/,7)  (2.5-21) 


where  ALAM  is  called  the  "energy  invariant";  it  is  defined  to  be  positive  for  positive  ions, 
negative  for  electrons.  For  an  isotropic  but  mono-energetic  particle  distribution,  this 
energy  invariant  remains  constant  as  the  particles  drift  along.  The  energy  invariant  ALAM 
is  the  analogy,  for  the  case  of  an  isotropic  plasma  distribution,  to  the  magnetic  moment  p. 
for  the  case  of  a  distribution  of  particles  that  all  have  90°  equatorial  pitch  angle  and  thus 
mirror  in  the  equatorial  plane.  In  (2,5-21),  VM  is  the  -2/3  power  of  the  flux  tube  volume 
(equation  (2.3-4)).  In  the  MSM,  VM  has  units  of  {R^nT)-'^  and  ALAM  has  units  of  eV 
(/?^nT)2/3.  For  a  proof  that  ALAM  is  the  appropriate  invariant  for  the  case  of  an  isotropic 
plasma,  see  Harel  et  al.  [1981]  or  H^o//  [1983].  The  interpolation  in  energy  required  to 
find  the  flux  at  energy  EilJJE)  is  done  by  straightforward  linear  interpolation  using  the 
log  of  the  energy  and  the  log  of  the  flux. 

25.4.2.  Calculation  of  the  Initial  Condition  Flux 


Let  jo{I,J,IE)  represent  the  initial-condition  flux  at  grid  point  (/,/)  and  energy 
channel  IE.  In  general,  we  calculate  this  flux  from  a  combination  of  FLXMAT  and  a 
previously  calculated  set  of  fluxes  for  the  same  time,  which  we  label  joidUJJE). 
Specifically,  we  write 

UU,  IE)  =  bojoid{U,  m —  (2.5-22) 

aoj\jKpi.U,  IE)]  +  boAioidUJ,  IE)] 


where  f(j)  =  1  if  j  >  0, 0  if  y  <  0. 

The  basic  logical  structure  of  the  MSM  allows  use  of  a  general  admixture  of 
FLXMAT  and  joid  in  setting  the  initial-condition  plasma  distribution.  However,  in  its 
current  configuration,  the  program  automatically  sets  Co  =  0  if  it  has  an  initial  distribution  to 
start  from,  and  bo  =  0  when  it  has  no  initial  distribution  to  start  from. 

2.5. 4. 3.  Conversion  of  Fluxes  to  Number  Invariant  rj. 

It  can  be  shown  that,  if  charged  particles  £xB-drift  and  gradient/curvature  drift 
adiabatically  and  losslessly,  then  the  number  of  panicles  per  unit  magnetic  flux  of  a  given 
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"species"  IE,  specifically  a  given  chemical  species  and  given  energy  invariant  ALAM,  is 
conserved  [Harel  et  al.,  1981;  Wolf,  1983].  Therefore,  we  define  r[{IE,IyJ)  to  be  the 
number  of  particles  per  unit  magnetic  flux  of  species  IE,  and  do  the  particle- trace 
calculations  in  terms  of  i]  rather  than  the  flux  j.  The  expression  relating  differential  flux 
and  T]  is  the  following: 

ieiectronil,  J,  IE)  =  7.392  X  lO'^^  X  VM{IJ)  X  r\{IEJJ) 

JUI.J.tE)=  ^2.5-24) 

A  [\'kmaxJE\  -  \1-minJE[ 


where 

XnunjEi^  0  for  the  lowest-energy  channel  for  a  given  chemical  species  (2.5-25a) 
^nunjE  =  otherwise.  (2.5-25/?) 


Similarly,  define 


“kmaxJE-  '2‘^IE  *  ^/£-1  (2.5-26fl) 

if  K  is  the  highest-energy  channel  for  either  electrons  or  ions,  and 

l^maxJE  =  otherwise.  (2.5-26/?) 

2.5.5.  The  Boundary-Condition  on  t]. 

It  is  easiest  to  discuss  the  boundary  condition  in  terms  of  the  number  invariant  T) 
rather  than  the  differential  flux  y.  We  define  ETABND(J,  IE,  UT)  to  be  the  bound^ 
value  of  T]  for  local-time  grid  point  J,  invariant  energy  channel  IE,  and  time  UT.  As  testing 
of  the  MSM  progressed,  we  found  that  the  system  performed  best  if  we  chose  ET ABND  to 
be  independent  of  J.  In  other  words,  we  are  assuming  that  T]  is  constant  everywhere  on  the 
boundary.  Bear  in  mind  that  the  boundary-condition  value  of  T]  influences  the  calculation 
only  for  boundary  points  on  which  the  particles  are  flowing  into  the  modeling  region  from 
the  boundary.  Assuming  T]  to  be  a  constant  on  the  boundary  is  a  conventional  assumption 
in  the  modeling  of  magnetospheric  convection.  The  effects  of  allowing  to  vary  on  the 
boundary  is  a  topic  of  current  research  interest,  but  the  research  results  are  not  yet 
sufficiently  clear  to  warrant  varying  T]  in  an  operational  model  like  the  MSM. 

We  choose  the  value  of  ETABND  by  the  requirement  that  it  correspond  to  FLXMAT 
at  a  point  in  the  equatorial  plane  that  is  at  local  midnight  and  13  Re  geocentric  distance.  In 
choosing  ETABND  to  match  FLXMAT ,  we  use  (2.5-23)  or  (2.5-24)  for  the  real-time  Kp, 
but  for  an  average  value  of  VM  for  a  13  Re  and  local  midnight,  namely  2.182.  We  chose 
to  match  ETABND  at  13  Re  .  rather  than  a  boundary  point,  because  the  major 
observational  statistical  statistical  studies,  like  that  of  Huang  and  Frank  [1986],  apply  to 
the  inner  plasma  sheet. 

This  procedure  for  choosing  ETABND  has  worked  well,  but  the  user  of  the  MSM 
should  be  aware  of  a  peculiar  characteristic  of  this  choice:  ETABND  remains  independent 
of  time  for  each  three-hour  Kp  interval,  then  shifts  abruptly  whenever  Kp  changes.  It 
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should  further  be  remarked  that  changes  in  the  boundary  value  of  Tj  do  not  usually  have 
dramatic  effects  on  fluxes  inside  synchronous  orbit,  because  of  our  use  of  a  reference  flux, 
as  discussed  in  the  following  section. 

2.5.6.  The  Reference  Flux. 

One  of  the  puzzles  of  present-day  magnetospheric  physics  is  the  apparent  violation  of 
the  adiabatic-convection  condition  in  the  inner  plasma  sheet  [Erickson  and  Wolf,  1980]. 
The  quantity  which  we  would  expect  to  be  approximate  adiabatic  invariant,  typically 
increases  strongly  with  increasing  distance  down  the  tail.  We  do  not  know  what  physical 
mechanisms  are  responsible  for  the  non-constancy  of  (For  a  discussion  of  this 

general  topic,  see,  e.g.,  Pontius  and  Wo// [1990].) 

For  the  MSM,  which  is  an  operational  model,  we  had  to  find  a  way  to  parametrize 
around  this  point  on  which  the  physics  is  not  understood.  If  we  assume  adiabatic 
convection  from  the  inner  plasma  sheet  (~13  Re),  we  obtain  fluxes  at  geosynchronous  orbit 
that  are  about  an  order  of  magnitude  larger  than  indicated  by  observations.  (See  Appendix 
g) 


We  take  care  of  the  problem  by  placing  a  ceiling  on  the  flux  values.  The  computed 
flux  values  at  all  grid  points  and  energy  channels  are  compared  with  a  reference 

flux  jr^r(lj),  E(I,J)).  If  the  computed  j(IEJJ)  >  jrej(r{IJ),  E{IJ)),  then  ;(/£, 7,7)  is  set 
equal  to  jrejiriU),  E{IJ)). 

To  compute  the  reference  flux,  we  first  compute  a  matrix  jmAr,E),  defined  by 

=  m&\{FLXMATiKp,r,E),  Kp=0,\ . 6}  (2.5-27) 

The  reference  flux  jrej(r(I,J),  E{I,J))  is  simply  set  equal  to  jmax{r{I,J),  E(I,J)),  for 
r(77)  ^  6.6.  For  r  <  6.6,  gradient/curvature  drift  prevents  direct  earthward  transport  of 
plasma-sheet  particles,  and  there  is  no  clear  need  to  invoke  a  "mystery  mechanism"  for 
particle  loss  inside  geosynchronous  orbit.  Therefore,  for  r  <  6.6,  we  set  yr^at  a  level  that 
would  correspond  to  adiabatic  compression  of  the  geosynchronous  flux.  Since  }ds/B  « 
in  a  dipole  field,  and  VM  «  r  *8/3,  we  write 

jrej{r,E)  =  j^Py^6.6,  (2-5-28) 

2.5.7.  The  Minimum  Flux. 

The  loss  algorithms  employed  by  the  MSM,  which  are  discussed  in  Section  2.6.3  and 
2.6.4,  all  set  the  loss  rate  equal  to  zero  for  a  given  species  when  the  flux  level  for  that 
species  falls  below  a  specified  minimum  level,  which  we  will  call  jmirSf.E).  This  minimum 
flux  is  defined  simply  by 

jmin{r,E)  =  min{FLXMAnKp,  r,  E),  Kp  =  0,  1,...6}  (2.5-29) 

2.5.8.  The  Threshold  Flux. 

An  additional  flux  matrix,  called  THRMAT,  is  of  the  same  general  form  as  FLXMAT 
and  is  calculated  along  with  it.  THRMAT  represents  the  threshold-flux  level  for  the  onset 
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of  strong  pitch-angle  scattering  for  electrons.  The  logic  for  the  calculation  of  THRMAT  is 
described  in  Section  2.6.3,  along  with  other  aspects  of  the  electron  loss  algorithm. 
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2.6.  Particle  Tracer  and  Loss  Algorithms 

2.6.1.  The  Equations  of  Adiabatic  Drift 

2.6.1. 1.  The  Differential  Equations. 

As  in  most  modem  efforts  at  modeling  magnetospheric  convection,  we  allow  particles 
to  drift  perpendicular  to  B  by  ExB  drift  and  gradient/curvature  drift.  We  further  assume 
that  particles  of  a  given  species  IE  (i.e.,  given  chemical  species  and  given  invariant-energy 
level),  are  isotropically  distributed  with  regard  to  pitch  angle.  In  this  case,  the  formula  for 
bounce-averaged  ExB  and  gradient/curvature  drift  is 

v/£  =  -  +  ALAM -  (2.6-1) 


where  ALAM  is  the  energy  invariant  for  an  isotropic  distribution  of  particles  in  adiabatic 
convection.  Our  use  of  ALAM  was  defined  in  Section  2.5.4. 1  and  is  discussed  in  detail 
by  Harel  et  al.  [1981]  and  Wo//[1983].  The  symbol  VM  represents  the  -2/3  power  of  the 
flux-tube  volume,  as  defined  in  Section  2.5.4. 1. 


The  assumption  that  the  particle  distribution  in  a  collisionless  plasma  remains 
isotropic,  combined  with  the  assumption  that  ALAM  is  conserved,  implies  that  the  particles 
undergo  frequent  wave-particle  interactions  that  result  in  pitch-angle  scattering,  but  not  in 
any  change  in  particle  energy.  It  is  the  best  simple  theoretical  description  of  the  Earth’s 
plasma  sheet,  which  is  genertdly  observed  to  be  highly  isotropic. 


The  electric  field  E  in  equation  (2.6-1)  in  general  includes  both  a  potential  field  and  an 
induction  field.  It  can  be  applied  at  any  point  on  a  field  line.  Out  in  the  magnetosphere,  the 
magnetic  Held  changes  substantially  in  time,  and  the  induction  electric  field  is  generally 
important.  In  the  ionosphere,  however,  the  changes  in  B  are  small  enough  to  be 
unimportant  from  our  point  of  view,  and  we  can  write  E  =  -  W,  where  V  is  the 
electrostatic  potential.  We  can  then  rewrite  equation  (2.6-1)  as 


where 


V/£ 


BxVVeff 

B^ 


(2.6-2a) 


V  +  ALAM  X  VM 


{2.6-2b) 


It  can  be  shown  that,  if  the  particles  conserve  ALAM  and  suffer  no  loss,  then  \\ie,  the 
number  of  particles  of  species  IE  per  unit  magnetic  flux,  is  conserved  along  a  drift  path. 
We  therefore  write 


P^tE  _ 
Dt 


(2.6-3) 


The  form  of  the  loss  term  will  be  discussed  in  Sections  2.6.3  and  2.6.4.  Equations  (2.6-2) 
and  (2.6-3)  are  the  basic  equations  that  the  particle  tracer  solves  numerically. 

The  particle-trace  procedure,  which  is  the  central  part  of  the  MSM,  is  a  useful  method 
of  follwing  electrons  in  the  geosynchronous-orbit  region,  for  particle  energies  below  about 
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100  keV.  Electrons  below  about  100  keV  are  strongly  affected  by  substorm  injections, 
whereas  most  variations  in  the  flux  above  -100  keV  are  not  obviously  substorm  related. 
There  is  an  analagous  effect  in  the  ions  near  geosynchronous  orbit.  The  top  panel  of 
Figure  2.5-4  suggests  that  fluxes  below  about  30  keV  vary  considerably  with  Kp,  but 
more  energetic  does  not.  Oxygen  fluxes  vary  significantly  with  Kp,  but  0+  tends  to 
exhibit  considerably  smaller  flux  levels  than  H+  for  energies  above  ~  50  keV. 

As  delivered,  the  MSM  is  set  not  to  perform  detailed  particle  traces  for  electrons  with 
geosynchronous  energies  above  about  100  keV  or  for  ions  with  geosynchronous  energies 
above  about  50  keV.  For  energies  above  those  levels,  the  fluxes  outputted  by  the  MSM  are 
computed  simply  by  scaling  the  appropriate  j^p  to  agree  with  the  observed 
geosynchronous  fluxes  for  the  previous  fifteen  minutes. 

The  decision  that  the  MSM  makes  on  whether  it  is  worthwhile  to  perform  detailed 
particle  traces  for  particles  of  a  given  invariant  energy  is  easy  to  adjust  within  the  MSM. 
However,  the  choice  that  is  presently  programmed  into  the  code  represents  our  present  best 
judgement  as  to  the  optimum  cutoff  point  for  the  tracing. 

2.6.12.  The  MSM  Coordinate  Grid. 


The  MSM's  ionospheric  grid  is  equally  spaced  in  local  time.  The  local-time  spacing 
is  7.5®,  which  corresponds  to  0.5  hr  of  local  time.  There  is  a  wraparound  of  3  grid  points, 
so  that  JMAX,  the  total  number  of  local-time  grid  points,  is  51. 


The  latitudinal  grid  spacing  is  non-uniform,  with  the  closest  spacing  in  the  latitude 
region  that  is  normally  the  auroral  zone  and  wider  spacing  in  the  polar  cap  and  the  low- 
latitude  region.  We  used  the  following  formula  for  the  latimdinal  grid  spacing: 


dl 


H 

dQ 

dd 

tan'^ 

\lpc-I 

•+■  tan'^ 

[dliruvc  tt  i 

dl 

max 

dl 

min 

A/ 

.  A/  J| 

(2.6-3) 


where  0  is  colatitude,  I  is  the  latitudinal  grid  index  (lowest  near  the  pole,  highest  near  the 
equator);  (dQldf)^,  {dQldI)^n^  A/,  I  pc,  and  are  adjustable  parameters  that  represent 
the  approximate  maximum  and  minimum  grid  spacing,  the  width  of  the  transition  region 
between  dense  and  undense  spacings,  and  the  /-locations  of  the  boundaries  of  the  dense- 
grid  region.  (Here  "pc"  means  "polar  cap",  "pp"  means  "plasmapause".)  Equation  (2.6-3) 
integrates  to  the  form 


0(/)  =  0(1)  +  (/-I)  M  + 

di  jmax 

.  A/  \dQ 

^  1 

where 


F{x)  =  xtan'Hjt)  - 


(2.6-4P) 


The  present  MSM  uses  the  following  numerical  values: 
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(d0/d/)^  =  2.0° 

(2.6-5a) 

(d0/d/)^,„  =  O.4° 

(2.6-5b) 

A/  =  2 

(2.6-5C) 

Ipc  =  6 

(2.6-5d) 

Ipp  =  55 

(2.6-5C) 

Figure  2.6-1  shows  a  plot  of  colatitude  vs.  I  for  the  coordinate  system  that  is  being 
delivered  with  the  MSM.  The  I  coordinate  runs  from  1  up  to  IMAX,  which  we  have 
chosen  to  be  62  in  the  present  case;  the  maximum  colatitude  is  49. 10°. 


In  the  program,  we  use  the  following  conventions  with  regard  to  grid-spacing: 


dl 


ALPHA  X  DLAM 


(2.6-6) 


where  [dQ/df\  is  in  radians,  and  the  constant  DLAM  is  equal  to  l/(/A/AX-l),  which  is  1/61 
in  the  present  case.  The  distance  between  between  local-time  grid  points 

BETAxDPSIxRI 


where  DPSf  =  2n/(JMAX-3)  and  RI  is  the  geocentric  distance  of  the  ionosphere,  which  is 
set  equal  to  6500  km.  It  follows  that  BET A  -  sin(6). 

The  RCM's  coordinate  grid  is  assumed  to  rotate  with  the  Earth.  At  midnight 
Universal  Time,  7  =  3  is  at  local  noon  in  the  ionosphere,  7  =  15  is  at  local  dusk,  7  =  27  is  at 
local  midnight,  and  7  =  39  is  at  local  dawn. 

The  colatitude  0  for  grid  point  {IJ)  is  called  COLAT(IJ)  in  the  program,  while  the 
local  time  of  the  same  grid  point  is  called  ALOCT(I,J).  The  general  structure  of  the 
program  allows  both  COLAT  and  ALOCT  to  vary  with  both  I  and  7.  However,  we  have 
not  taken  advantage  of  that  flexibility  in  the  present  version  of  the  MSM:  at  present, 
COLAT  really  depends  only  on  I  and  AZ79CT  depends  only  on  7. 

The  present  version  of  the  MSM  does  not  include  the  effects  of  dipole  tilt,  as 
discussed  in  Section  2.3.  The  code  implicitly  assumes  that  the  Earth's  magnetic  dipole  is 
aligned  with  the  rotation  axis.  Also,  the  magnetic-field  strength  at  ionospheric  height  is 
assumed  to  be  purely  dipolar.  For  a  planet  whose  magnetic  field  is  purely  dipolar  and 
whose  dipole  axis  is  aligned  with  the  rotation  axis,  geographic  coordinates,  geomagnetic- 
dipole  coordinates,  corrected  geomagnetic  coordinates,  and  invariant  latitude  all  become 
coincident.  When  applying  model  predictions  to  the  real  Earth,  it  is  best  to  interpret  the 
model's  COLAT  as  either  invariant  colatitude  or  corrected  geomagnetic  colatitude  and  the 
model's  ALOCT  as  magnetic  local  time. 
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2.6.13.  Equation  of  Motion  in  Terms  of  the  MSM  Grid. 

The  drift  equation  (2.6-2a),  written  in  terms  of  motion  of  a  particle  in  the 
system,  takes  the  following  form: 

^  _  _ 1 _ 

ALPHAxBETAxDLAMxDPSlxBIRxRl^  I  1/ 


1 

(3V 

ALPHAxBETAxDLAMxDPSlxBlRxRl^ 

1  a/  j 

MS  M's  grid 


(2.6-7a) 


(2.6-7b) 


where  BIR  is  the  strength  of  the  radial  component  of  the  magnetic  field  at  ionospheric 
altitude  (taken  to  be  positive). 

2.6.2.  Operation  of  the  Particle  Tracer. 

2. 6. 2.1.  Walking  Test  Particles  Backwards  in  Time. 

Suppose  that  we  know  the  flux  levels  at  all  grid  points  (IJ)  and  for  all  invariant- 
energy  levels  /£,  for  time  r,  and  we  wish  to  step  along  to  time  t  +  At.  (The  basic  time  step 
At  of  the  model  is  set  in  the  code  to  be  15  minutes;  our  choice  of  At  was  arbitrary,  and 
could  easily  be  changed.) 

To  calculate  the  plasma  distribution  at  time  t  to  calculate  the  known  distribution  at  t  + 
At,  the  particle-trace  program  utilizes  the  electric-  and  magnetic-field  configurations,  which 
were  computed  for  all  times  through  the  run  before  the  particle  tracer  was  first  called. 
Using  this  known  fields,  the  panicle  tracer  starts  a  test  particle  of  each  species  IE  from 
each  grid  point  (IJ)  and  walks  backward  in  time  from  r  +  At  to  t.  The  objective  of  the 
backward  walk  is  to  determine  the  location  Uuji)  of  the  test  particle  at  time  t. 

The  backward  walk  is  accomplished  numerically  by  means  of  a  fourth-order  Runge- 
Kutta  procedure,  with  a  fifth-order  correction.  The  Runge-Kutta  time  step  is  adjusted 
automatically  by  a  standard  error-check  algorithm.  When  it  needs  to  calculate  the  velocities 
(equation  (2.6-7))  at  an  arbitrary  point  {i’,j')  and  time  t\  the  program  calculates  ALPHA, 
BETA,  and  BIR  by  straightforward  linear  interpolation  amongst  the  neighboring  grid 
points.  To  calculate  dV^dl  and  dVgHdJ,  the  program  linearly  interpolates  in  time  and 
uses  central  differences  to  calculate  oVg^dl  and  dVgj^dJ  at  the  four  grid  points  nearest 
(i',y')-  It  then  calculates  the  derivatives  at  (/',;')  by  straightforward  2d  linear 
interpolation. 

2. 6. 2. 2.  The  Loss  Equation  and  its  Integration. 

The  algorithm  described  here  divides  the  loss  into  two  types,  one  that  involves 
strong  pitch-angle  scattering  above  a  threshold,  another  that  involves  a  lower  level  of  loss 
due  to  some  kind  of  background  noise. 

We  assume  that  loss  is  governed  by  the  following  differential  equation: 

Dr{  ^  Max[0,(ii--n,;,rg5ii)]]  ^2  6-8) 

Ot  L  "^weak  '  “^strong  J 
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where  Max(j:,}')  is  the  laiger  of  x  and  _y,  (XweaJt)'^  is  a  background-noise-related  loss  rate 
that  is  far  below  the  strong-pitch-angle  scattering  rate  (discussed  in  Section  2.6.3. 1),  and 
i'^strong)'^  is  a  rate  that  is  close  to  the  strong-pitch-angle-scattenng  limit  (discussed  in 
Section  2. 6.3. 2).  The  parameter  r]:/,resh  is  a  threshold  density  for  the  onset  of  strong  pitch- 
angle  scattering  (discussed  in  Section  2.6.3. 1).  It  is  additionally  assumed  that  there  is  a 
threshold  for  weak  pitch-angle  scattering:  it  is  assumed  that  weak  scattering  drives  the 
invariant  density  down  to  a  low  level  r\eq,  which  corresponds  to  jmin  ,  as  defined  in 
equation  (2.5-29).  Imposition  of  the  loss  cutoff  Tj^^,  prevents  the  invariant  density  from 
dropping  to  unphysically  low  values. 

The  loss  calculation  is  carried  out  after  the  traceback  for  a  giv'^n  test  particle  from 
r-t-Ar  to  time  t  has  been  completed.  Specifically,  we  assume  that  the  program  remembers, 
for  each  step  m  of  the  traceback,  the  values  of  the  weak  loss  rate  (twea*  ')  and  the  strong 
loss  rate  (tsirong  *)  at  the  beginning/end  of  each  time  step.  The  program  also  remembers  the 
density  levels  x\eq  and  T[thresh  at  the  interval  end  time  r-t-A/. 

Consider  now  the  decrease  in  q  during  the  time  step  from  tm-i  to  tm,  which  is 
between  mark  time  t  and  the  next  mark  time  t+d^t.  We  first  check  to  see  which  of  the  two 
terms  on  the  right  side  of  (2.6-8)  is  to  be  taken  as  larger  for  this  step.  We  do  this  using 

Max[0,('n(m-l)--ng^(/))] 

for  the  first  term  and 

Max[0  ,(q  (m- 1  yx\threshif) )] 


for  the  second;  here  we  take 


t/(m-.5)  =  ^[t,(m-l)  t,(m)] 

where  i  is  either  "weak"  or  "strong".  If  both  of  these  terms  are  zero,  then,  of  course,  T| 
remains  constant  from  to  tm,-  If  at  least  one  is  positive,  then  the  differential  equation 
(2.6-8)  takes  the  form 


^  (2.6-9) 

Dt 

The  parameters  "n,  and  Xm  are  taken  to  be  either  ^^(^(t-i-Ar)  and  iweakim-.5),  or 
^thresh(t-^^0  and  Xstrong{fn-.5),  whichever  is  appropriate.  In  either  case,  we  take  r\,_m  and 
Xm  to  be  constants  through  the  time  step,  which  means  that  (2.6-9)  can  be  immediately 
integrated  to  give 


q(m)-Ti,(m)  =  [q(m-l)  -  q,(m)]  exp|  -  (2.6-10) 

We  walk  forwards  in  time  through  the  interval  from  /  to  t+At  by  repeated  use  of  (3). 

2.6.3.  Electron  Loss 

2.6.3. ] .  The  Electron  Strong-Pitch-Angle-Scattering  Rate 
The  strong-pitch-angle-scattering  loss  rate  is  defined  by 
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yuii 


(2.6-11) 


where  v  is  the  velocity  of  the  particle’s  gyro-  and  bounce  motion,  and  is  the  magnetic- 
field  strength  at  the  ionosphere.  Assuming  loss  at  the  rate  given  by  (2.6-1 1)  corresponds 
to  assuming  that  the  pitch  angle  is  everywhere  isotropic,  so  Aat  the  loss  cone  is  completely 
full.  Rewriting  (2.  1)  in  the  terminology  and  units  of  the  MSM  gives 


- 1 -  =  (4.45  X  (2.6-12) 

yuii[x(t)] 

If  X,  VM,  and  Bir  are  in  our  usual  program  units  [  eV  (Rf/nT)^/^,  (/?£/nT)-2/3,  and  nT, 
respectively]  and  m  is  in  kg,  then  this  loss  rate  comes  out  in  sec  ^  The  base  loss  rate  is 
given  by  this  full  strong-pitch-angle-scattering  rate,  multiplied  by  an  additional  factor  that 
we  estimate  on  the  basis  of  results  obtained  by  Schumaker  et  al.  [  1989].  Specifically,  we 
set 


where 


T-i—  =  (2.6-13) 

^strong  yfuU 

fo  =  0.333  for  ^:p<  1.5  (2.6-14fl) 

fo  =  0.333  i2.5-Kp)  +  0.667  (Kp-1.5)  for  1.5  <Kp<  2.5  (2.6-14f?) 

fo  =  0.667  for  Kp>  2.5  (2.6- 14c) 


2.6.3. 2.  jihreshoid  for  Electron  Strong  Pitch-Angle  Scattering. 

The  classic  paper  of  Kennel  and  Petschek  [1966]  predicted  that  there  would  be  a 
threshold  flux  for  the  onset  of  electron  strong  pitch-angle  scattering  by  whistlers,  and  that 
that  threshold  would  be  approximately  proportional  to  L^-.  On  the  other  hand,  the  paper 
by  Baker  et  al.  [1979]  found  observational  evidence  for  such  a  threshold,  but  set  it  at  ~  5  x 
10^  cm-2  S'*  sr*  for  £  >  40  keV,  which  is  about  10  times  the  level  originally  suggested  by 
Kennel  and  Petschek.  We  thus  write 


•/ifPrfe-«A(>40keV)  =  (5  x  10  ^  cm-2  s'* 


sr*) 


r  J 


(2.6-15) 


Equation  (2.6-15)  can  be  related  to  differential  flux  for  a  kappa-type  distribution  function, 
using  equation  (2.5-2).  Integrating  the  resulting  expression  over  energy  from  40  keV  to 
infinity  gives 


y(>40  keV)  =  (1.68  X  10*  cm-2  s'*  ster*)  £(k) 


1/2 

N  1 

1  1  40keVl 
L  Eo 

ll  keVJ 

.1  cm-3. 

[l+40keV]’' 

kEo  J 


(2.6-16) 
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We  thus  define  an  energy  ratio  Er  as  follows: 


y(>40  kev) 


E_ 

Eo 


/n-4QkeV 

1  ^Eq 


(i 

\  kEJ  ' 


+4QkgV\ 
E 


I 


(2.6-17) 


Note  that,  if y(£)  is  to  come  out  in  cm'^  s*^  eV'^  sr^  as  is  our  usual  convention,  then,  in 
the  initial  multiplying  factor  \/Eo  on  the  right  side  of  (14),  Eo  should  be  in  eV. 


For  r  =  3, 4,  6.6,  and  13,  we  can  use  the  values  of  k  and  Eo  that  were  adopted  in  the 
calculation  of  FLXMAT,  as  described  in  Section  2.5.2.  They  are  summarized  in  the 
following  table: 


r 

Kp 

K 

£«(keV) 

5 

all 

2.8 

23.5 

4 

all 

2.0 

16 

6.6 

(i 

6 

21 

6.6 

1 

6 

4.8324 

6.6 

2 

6 

4.8140 

6.6 

3 

6 

4.7978 

6.6 

4 

6 

4.8160 

6.6 

5 

6 

4.8943 

6.6 

^6 

6 

5.0752 

13 

aU 

6 

.2885(/rp-l)-i-.15625(5-/(:p) 

To  calculate  the  threshold  flux,  then,  we  use 


jihreshir,  E)  =  (5  X  10  cm-2  S'l 


sri)  £/j(r,£) 


(2.6-18) 


Equation  (2.6-18)  is  designed  to  give  the  threshold  flux  for  energies  above  about  40 
keV.  Below  20  keV,  the  results  of  Schumaker  et  al.  [1989]  indicate  that  the  scattering  rate 
is  Close  to  the  strong  limit,  as  indicated  in  (2.6-13)  and  (2.6-14),  but  they  do  not  suggest 
evidence  of  a  clear  threshold.  However,  to  prevent  fluxes  from  dropping  to  unphysical 
low  values,  we  continue  in  the  MSM  to  use  a  threshold  value  of  r\eq  that  represents  a 
quiet-time  minimum  flux  level. 

To  span  our  required  energy  range,  I  propose  that  we  use  a  mixed  algorithm.  For 
energies  below  20  keV  or  below  40  keV,  we  make  the  obvious  choices: 


j thresh  =jeq  for  £  <  20  keV 


(2.6-19a) 


jthresh  =  value  from  (2.6- 1 8)  for  £  >  40  keV  (2.6- 1 9b) 


Between  20  and  40  keV,  we  do  a  log-linear  interpolation  in  energy: 


^O^threshiE))  — 


(40-£)  lo^'thresh(20))  +  (£-20)  lo^threshW) 
20 


(2.6-20) 
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for20keV<£<40keV. 

Near  the  beginning  of  each  run,  the  MSM  computes  a  matrix  of  values  of 
^ogioijihresh),  in  analogy  with  the  matrix  FLXMAT.  This  provides  values  of  logioOWcj/i) 
for  Kp=0,  1,  2,  3,  4,  5,  and  6,  for  r  =  3,  4,  6.6,  and  13,  and  for  appropriate 
logarithmically  spaced  energy  levels.  As  usual,  we  can  convert  jthresh  to  T[ihresh  using 
equation  (2.5-23). 


2.63.3.  Weak  Loss  Rate 

Calculation  of  the  weak  loss  rate  depends  on  whether  we  are  inside  or  outside  of  the 
plasmapause.  Although  we  don't  keep  track  of  the  plasmapause  in  the  model,  we  can 
estimate  its  equatorial  geocentric  distance  as  follows: 


tpp  - 


(2.6-21) 


inside  the  plasmapause.  If  r  <  r^p,  then  we  directly  use  Lyons'  classic 
calculation  of  the  loss  lifetimes  of  electrons  due  to  pitch-angle  scattering  from 
plasmaspheric  hiss  Lyons  et  al.  [1972].  We  found  log-parabolic  analytic  approximations 
to  the  published  curves.  The  results  for  20, 50, 200,  and  500  keV  are  as  follows: 

logio(Tw^aik(20keV))  =  (  4.9365  +  15.8621)  -  6.06897  r  +0351724  r2 

(2.6-22a) 

logio(tw^ai(50keV))  =  (  4.9365  +  133172)  -  538621  r +0351724 

{2.6-22b) 

logio(TwMjt(200  keV))  =  (  4.9365  +  1.15517)  -  0.42069  r  +0.0275862 

(2.6-22C) 


logio(TH,«a*(500  keV))  =  (  4.9365  +  3.89655)  -  1.75862  r  +0206S97 

(2.6-22d) 

For  energies  below  50  keV,  use  a  linear  interpolation  between  (2.6-22a)  and  (2.6-22^?);  for 
energies  between  50  and  2(X)  keV,  use  linear  interpolation  between  (2.6-22fc)  and  (2.6- 
22c),  etc.  Of  course,  we  shouldn't  really  calculate  fluxes  for  energies  above  about  300 
keV,  since  we  aren't  even  using  relativistic  drift  formulas.  Figure  2.6.2  shows  a  plot  of 
lifetimes  computed  from  these  algorithms. 

Estimation  of  outside  the  nlasmapau.se.  We  have  very  little  solid  information 
concerning  x^eak  outside  the  plasmapause,  but  we  have  to  make  a  reasonable  estimate.  It  is 
known  that  a  substantial  fraction  of  the  ~  50  keV  electrons  are  lost  as  they  traverse  the  day 
side  of  the  magnetosphere,  which  takes  ~  an  hour.  We  therefore  assume  that,  for 
synchronous  particles  with  50-300  keV,  the  peak  dayside  loss  lifetime  is  ~  0.5  hr.  Based 
on  the  plots  in  the  paper  by  Frank  et  al.  [1964],  we  estimate  that  the  peak  dayside 
precipitation  rate  is  -  10  times  the  nightside  rate.  Therefore,  we  assume  that 
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('Cw<.a/fc(6.6))-^  =Min(— £ - l.ll  ( — 1 — +(-^)cos(<l>) 

^20keV  /l3600s/L  111/ 


for  £  >  20  keV 


(2.6-23a) 


(Wflt(6.6))'*  =0  for£<20keV  (2.6-236) 

There  is  no  need  for  this  extra  term  for  electrons  below  20  keV,  because  the  strong-pitch- 
angle-scattering  threshold  is  low  there  anyway. 

For  geocentric  distances  between  rop  and  6.6,  we  use  a  linear  interpolation  between 
Lyons'  results  for  the  plasmapause  and  (f.6-23),  as  follows: 

(W(r)^‘  (2-6-240) 

Of  course,  if  rpp  >  6.6,  we  shouldn't  use  (2.6-24fl).  For  r  >  6.6,  assume  that 

{W(r))-^  =  (tl!°i*(6.6))'*  (2.6-246) 


2.6.4.  Ion  Loss 

2.6.4.O.  Introductory  Comments. 

Loss  of  magnetospheric  ions  takes  place  by  two  physically  different  mechanisms: 
precipitation  and  charge  exchange.  Because  charge  exchange  is  more  significant  overall 
and  can  be  estimated  much  more  reliably,  we  have  simply  neglected  ion  loss  by 
precipitation.  In  applying  equation  (2.6-8),  charge  exchange  is  treated  as  "weak  loss". 
The  strong-loss  rate  is  set  equal  to  zero  for  ions. 

The  main  technical  discussion  of  charge-exchange  loss,  Section  2.6.4. 1,  was  written 
by  Dr.  James  Bishop  (Department  of  Atmospheric  and  Oceanic  Science,  University  of 
Michigan,  Ann  Arbor,  MI  48109),  who  performed  the  calculations.  Technical  questions 
about  the  calculations  should  be  directed  to  him.  His  phone  number  is  (313)936-0516. 

Section  2.6.4.2  describes  the  tabular  form  in  which  the  charge-exchange  calculations 
are  included  in  the  MSM. 

2. 6.4.1.  Charge  Exchange  Theory. 

We  consider  a  stationary  situation,  in  which  ion  energy  (speed)  remains  constant 
along  the  bounce  path.  The  degree  of  erosion  is  determined  solely  by  the  amount  of 
geocoronal  atomic  hydrogen  "experienced"  by  the  ring  current  ions  -  i.e.,  the  column 
depth  along  the  helical  path  executed  by  the  mirroring  ions.  Note  that  charge  exchange 
((2E)  collisions  act  strictly  as  a  loss  mechanism.  (We  will  neglect  the  thermal  atom  speeds.) 

We  will  focus  on  H+  and  Q*’  ions  with  energies  1-200  keV. 

First  consider  the  ideal  case  where  the  initial  equatorial  kinetic  distribution  function 
(KDF)  for  ring  current  (RC)  ions  of  species  i  is  isotropic  and  mono-energetic: 
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Mv,a,Qm  -  v)  (2.6-25) 

4nv^ 

where  niis.t)  is  the  density  of  species  i  at  location  s  (arclength  along  the  field  line  from  the 
equator,  sometimes  we  will  use  magnetic  latitude  A)  and  at  time  t,  andfi{v,a,s,t)  is  the 
underlying  KDF  at  speed  v  and  pitch  angle  a.  Further,  we  assume  that  initially 
fi(v,a,s,0)  is  equivalent  to  /i(t;,a,0,0)  in  the  sense  of  Liouville’s  equation  -  i.e.,  /i  is 
constant  dong  the  bounce  path  -  and  that  no  "injection"  occurs  for  r  >  0.  Then  we  can 
write  for  the  subsequent  evolution 

/i(v,a,0/)  = /(vja,0j0)  e-''i(“.0  (2.6-26) 


where  ti(a,r)  is  the  "depth"  of  geocoronal  H  traversed.  ti(a,r)  is  a  trajectory-dependent 
quantity,  conveyed  in  the  notation  by  the  dependence  on  pitch  angle  a: 

^1(0,0  =  J dlniiOiiv)=  dt'v(t)  «h(KO)  ai(u(t'))=  v  a-^iv)  |  dt'nii  (r(r')) 

(2.6-27) 

where  dl  is  an  element  of  path  length  along  the  ion  trajectory.  We  have  used  the  constancy 
of  V  and  the  fact  that  Oi  (the  charge  exchange  cross  section  for  ion  species  i)  is  a  function 
of  speed  alone.  The  equatorial  density  is 

rtiiOf)  =  I  d^  t>/i(t7,a,0/) 

Jdtpj  da  sin  a  I  dv'v^  ^i(Ofi)  e-Ti(a.o 

0  Jo  Jo 


=  ni(0,0)  I  da  sin  a 


(2.6-28) 


Several  comments  are  required: 

(i)  In  equation  (2.6-27),  we  have  used  the  assumption  of  Liouville  equivalence  along  the 
bounce  path  at  time  t  =  0  and  time  reversal  symmetry:  An  assemblage  of  ions  with  speeds 
in  the  range  (v,  v+dv)  and  pitch  angles  in  the  range  (a,a-t-da)  passing  the  equator  at  time 
t  =  0  will  all  be  located  at  time  (  at  a  position  s'  with  speeds  (y',  v'+dv')  and  pitch  angles 
(a',a'+da')  in  the  absence  of  CE  erosion.  An  assemblage  initially  at  s'  with  speeds  {v\ 
v’+dv')  and  pitch  angles  (7i-a',  7t-(a+da'))  will  at  time  f  be  located  at  the  equator.  Each 
assemblage  traverses  the  same  helical  column  of  geocoronal  H,  so  each  will  experience  the 
same  erosion. 


(ii)  The  quarter-bounce  time  (the  time  for  an  ion  of  speed  v'  and  pitch  angle  a  to  go  from 
equator  to  mirror  point)  is 
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for  a  dipole  magnetic  field;  here  B  is  the  magnetic  field  strength,  rgq  is  the  equatorial  radius 
of  the  field  line,  is  the  arclength  along  the  field  line  to  the  mirror  point,  and  Ab  is  the 
magnetic  latitude  of  the  mirror  point  (the  subscript  B  denotes  mirror  point  quantities).  The 
integral  is  speed  and  species  independent.  The  helical  quarter-bounce  column  density  of 
geocoronal  H  is 


The  total  column  of  geocoronal  H  traversed  by  an  ion  of  speed  t?  in  a  duration  t  can 
be  written 


1'*' 


'IhCKO)  =mT%{r, 


=  (m+l) 


(' 

eqitXeq)  +  I 

Jo 

Jr’O) 


ds'  ^ —  ,  m  even 

vnis')/v 


ds'  —  ,  m  odd 


(2.6-31) 


where  m  is  the  number  of  full  quarter-bounces  executed  in  the  duration  t  and  s{t)  denotes 
the  location  along  the  field  line  at  the  time  t.  If  t »  tb,  then  to  good  approximation  we  can 
drop  the  "remainder"  terms  and  work  entirely  in  terms  of  quarter-bounce  quantities.  In 
particular,  we  can  evaluate  t^XB(''eq.oteq)  and  NB(''eq.cteq)  at  the  outset.  The  evolution  is 
then  readily  followed  using 


ti(a,0  =  m  Oi  (V)  N B(''eq,aeq)  (2.6-32) 

where  m  =  INT(//'tB);  note  that  m  is  a  function  of  v,  Oeq,  and  r^q. 

(iii)  In  practice,  pitch  angles  are  limited  by  the  loss  cone  atoss-  ions  with  pitch  angles 
a  <  OLioss  have  mirror  points  so  deep  in  the  atmosphere  that  they  are  effectively  removed 
on  the  first  quarter-bounce.  Thus,  when  the  a-integral  is  evaluated,  the  summation 
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algorithm  will  generate  quadrature  points  dependent  on  aiosSy  which  in  turn  is  dependent 
onreq. 

(iv)  The  equatorial  density  should  track  with  the  flux  tube  content  as  long  as  xb  «  Xeff, 
where  the  effective  decay  time  for  equatorial  pitch  angle  <Xeq  is 

(Teffr^=^  ^  =  ^^^^  =  vai(v)  rniirit))  (2.6-33) 

fi  dt  dt 

where  r{t)  is  a  function  of  (and  r^a).  Recognizing  vxb  to  be  the  helical  quarter-bounce 
path  length,  a  mean  /i//-value  is  N  B/tnB,  so  a  mean  x^g^-value  for  the  bounce  motion  is 
a,(t;)NB/XB,.  Thus,  as  long  as  1  »  ai(t;)NBm  we  can  use  the  equatorial  density  as  a 
measure  of  flux  tube  content.  For  typical  geocoronal  conditions  <  10  *^  cm^,  so  we 
require  O  «  lO"^^  cm^,  which  is  generally  the  case.  (Note  that  another  problem  arises  if 
the  durations  of  evaluation  t  are  comparable  to  Xb,  since  then  we  cannot  neglect  the 
remainder  terms  in  (2.6-31).) 

Now  we  consider  the  case  where  /,•  has  a  binned  speed  distribution,  but  is  still 
initially  isotropic; 


/i(0 ,  a,  0,0)  =  !i<M  H(Vi.  Vi^i)  hi  (2.6-34) 

itl 

where  i;,-,  /  =  1,  2,...,  N  are  the  speeds  defining  the  N-\  bins,  H{Vi,  Vi+i)  =  1  for 
Vi  <  Vi,+\  and  0  otherwise,  and  the  weights  hi  for  each  bin  are  normalized  so  that 


Thus  we  consider  a  discrete,  finite  set  of  speeds  representative  of  each  bin;  CjJ-  1,  2,..., 
N-\.  The  equatorial  density  in  this  case  is 


/ii(0,r) 


where 


N  /”” 

^vfiiv,  a,  0,  0  =  n/(0,0)  X  I  ' 

/=]  Jo 


da  sin  a  (2.6-35) 


Xij  ~  fn  Oi{Cj)^g{rgq,  aeq) 


(2.6-36) 


Please  refer  to  Figures  2.6-3  to  2.6-5  for  sample  numerical  results  of  the  calculation. 

With  regard  to  accuracy,  the  primary  source  of  error  in  the  basic  charge-exchange 
calculation  probably  lies  in  the  use  of  a  very  simple  model  exosphere  and  in  using  a  crude 
approximation  to  the  solar-cycle  dependence  of  that  atmosphere.  For  L  <  4,  where  charge 
exchange  is  fast,  we  hope  that  the  rates  are  accurate  to  a  factor  of  two.  For  large  L,  the 
assumption  of  a  ^pole  field  in  the  calculation  decreases  accuracy  somewhat  further. 

2.6.42.  Incorporation  of  Charge -Exchange  Calculations  in  the  MSM. 
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In  the  basic  charge-exchange  calculations,  each  ion  is  assumed  to  drift  on  its  original 
L-shcU.  However,  the  ion  pitch-angle  distribution  changes  in  time,  due  to  the  rapid  erosion 
of  particles  with  small  pitch  angles,  which  experience  higher  average  neutral  densities  than 
particles  with  pitch  angles  near  90°. 

It  is  clear,  however,  from  Figure  2.6-5  that,  after  a  relatively  brief  interval,  the  overall 
density  for  given  E  and  L  decays  approximately  exponentially.  Since  the  implementation 
of  the  charge-exchange  calculation  in  the  MSM  has  to  execute  extremely  fast,  we  chose  to 
take  advantage  of  this  approximate  exponentiality,  to  avoid  the  necessity  for  having  the 
MSM  keep  track  of  the  time  evolution  of  of  the  pitch-angle  distribution  on  each  flux  tube. 
Dr.  Bishop  provided  us  with  a  table  of  loss  lifetimes  for  and  0+  for  various  energies 
and  L-values,  and  also  for  both  high  and  low  solar  activity.  The  MSM  reads  this  table 
from  disk  near  the  beginning  of  each  run.  Charge-exchange  rates  are  computed,  time- step 
by  time-step,  by  interpolation  of  the  table  in  energy  and  L-value. 


2.6.5.  Calculation  of  the  Fluxes  of  Precipitating  Particles 

2.65.0.  Introductory  Comments. 

Ion  precipitation  in  the  MSM  is  treated  simply  in  terms  of  the  statistical  model  of 
Hardy  et  al.  [1989],  which  is  the  most  reliable  model  available  at  present.  A  subroutine 
embodying  tat  model  was  kindly  provided  by  D.  A.  Hardy  of  the  USAF  Geophysics 
Laboratory. 

Within  the  region  where  we  carry  out  detailed  particle  traces,  we  compute  the  energy 
flux  in  auroral-electron  precipitation  in  terms  of  loss  from  our  model  plasma  sheet. 
Electrons  are  assumed  to  precipitate  at  the  rate  calculated  according  to  the  algorithms 
described  in  Section  2.6.3,  with  no  acceleration  by  field-aligned  potential  drops.  The 
procedure  for  calculating  the  energy  flux  and  average  energy  of  the  precipitating  electrons 
in  described  in  Section  2.6.5. 1. 

Poleward  of  our  main  modeling  region,  we  use  a  procedure  that  is  based  indirectly  on 
the  statistical  model  of  Hardy  et  al.  [1985].  It  is  described  in  Section  2.6.5. 2. 

2.65.1.  Calculation  of  Precipitating  Electron  Fluxes  within  our  Main  Modeling 
Region. 

We  compute  the  auroral  precipitation  from  the  invariant  densities  that  are  computed 
by  the  main  particle-trace  procedure,  subjected  to  loss  as  discussed  in  Section  2.6.3. 

For  each  (lEJJ),  it  is  useful  to  define  a  net  rate  (which  is  in  sec  '  and  is  the  inverse 
of  the  mean  precipitation  lifetime  of  the  particles  in  question): 

RATE  =  Max[0,(T)-Ti,;„e,/,)] 

fl  1  '^weak  '^strong 

The  number  of  particles  of  species  IE  precipitating  from  a  flux  tube  of  unit  magnetic  flux 
into  a  unit  area  in  the  ionosphere  (one  hemisphere)  is  given  by 


(2.6-37) 
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RATE*l\*Bir 

2xl0‘^ 


(2.6-38) 


where  the  factor  of  2  comes  from  the  fact  that  C>^r  is  a  one-hemisphere  rate,  and  the  factor 
of  10^3  comes  from  converting  from  nanoteslas  to  teslas  and  from  m-2  to  cm‘2.  Oyy 
comes  out  in  particles/cm^/s.  The  precipitating  energy  flux  per  hemisphere  in  species  IE  is 
given  by 


<D£  =  {0.'&xl0-^)RATE*T]*Bir*\'k\*VM  (2.6-39) 


where  the  constant  has  been  multiplied  by  1.6x10  *2  to  convert  from  eV  to  ergs.  The  flux 
Of  is  in  ergs/cm2/s  for  one  hemisphere  and  one  species  IE.  The  total  hemispheric  energy 
flux  for  one  chemical  species  is  given  by 


FLXSUMiUJK)  =  X  ^eUEJJ)  (2.6-40) 

IE 

where  the  sum  over  IE  includes  all  IE  -values  that  correspond  to  chemical  species  IK.  The 
corresponding  total  number  flux  is  given  by 

DENSUM(UJK)  =  X  ^nUEJJ)  (2.6-41) 

IE 


The  average  energy  is  given  by 


EAVGiUm  = 


FLXSUMiUJK) 

DENSUMiUJK) 


(2.6-42) 


2.65.2.  Electron  Precipitation  Poleward  of  our  Main  Modeling  Region. 


Our  formula  for  electron  energy  flux  is  expressed  in  terms  of  P"(I,J),  the  power  in 
precipitating  electrons  within  the  grid  space  (7,7),  in  units  of  gigawatts.  The  relation 
between  P"(U)  and  E/iux,  the  precipitating  energy  flux  in  ergs/cmvs,  is 

F’iU)  =  10-<*  Eflux  ail)  p(/)  AX  Av  /?/  (2.6-43) 

P"(iJ)  is  assumed  to  be  a  cubic  function  of  i,  the  floating-point  version  of  the  usual 
ladtudinal  grid  index,  and  specifically  to  have  the  form 


F'iiJ)  =  a  [(/-//)  +  b  U-ii)^  +  c  ii-ii)^  ]  (2.6-44a) 


for  i  >  ii  ,and 


F'iU)  =  0  il.e-AAb) 

for  i  <  ii .  Here  U  is  the  /-location  of  the  intersection  of  ellipse  #1  with  grid  line  J.  Note 
that  P"  is  forced  to  go  to  zero  at  that  ellipse.  The  three  coefficients  a,  b,  c  are  determined 
by  the  following  requirements: 
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(i)  the  cubic  formula  must  agree  with  the  results  of  the  low-latitude  computation  for  the 
first  two  grid  points  I  that  occur  within  the  main  modeling  region,  for  local-time  grid  line 
J. 

(ii)  it  must  give  a  designated  amount  of  power  (fj)  per  unit  J  poleward  of  the  main 
modeling  boundary.  For  diat  power,  we  use  the  algorithm: 


Pt 


JmIkiiN 


di  F'ilJ)  = 


F 

id(Power) 

Li-pJ 

\  dl 

(2.6-45) 


in  main  modeling  region 


where 


F  =  0.5  +  0.3  sin((l>) 


(2.6-46) 


and  ^  is  our  local-time  angle,  (0  at  noon,  nil  at  dusk....).  Formula  (2.6-46)  was 
originally  derived  by  comparing  the  statistical  electron-precipitation  patterns  of  Hardy  et  al. 
[1985]  with  the  Birkeland-current  patterns  of  lijima  and  Potemra  [1978],  associating  the 
boundary  of  our  main  modeling  region  with  the  equatorward  edge  of  the  region-one 
currents.  F.  J.  Rich  of  the  USAF  Geophysics  Lab  is  carrying  out  a  much  more  detailed 
study  of  the  question  of  where  the  region-one  current  boundary  falls  with  respect  to  the 
precipitation  pattern;  his  preliminary  results  are  reasonably  consistent  with  (2.6-46). 

Determination  of  ir.  The  size  and  shape  of  ellipse  #1,  which  is  the  poleward  edge  of 
the  electric-field  reversal  region  of  the  ionosphere,  is  defined  by  the  parameters  /1(1)  R(l), 
DX(l),  and  Dy(l).  From  the  equation  of  the  ellipse,  we  calculate  the  colatitude  at  which 
that  ellipse  crosses  local-time  grid  line  /.  We  then  find  the  floating-point  /-value  at  which 
the  ellipse  intersects  local-time  grid  line  J  by  simple  linear  interpolation  on  the  COLAT 
matrix. 

Determination  of  a.  b.  and  c.  Let  Pi"  and  Pi'  represent  P”(IMIN+\,J)  and 
P"(IMIN-¥2J),  respectively,  where  IMIN  is  defined  to  be  the  largest  /-value  that  lies 
poleward  of  our  main  modeling  region.  Also  let  iz  =  IMIN  -  i/. ,  /’i  =  IMIN  +  1  -  //, , 
i2  =  IMIN +2  -  ii.  The  formulas  then  are 


ii  -  Pz"  +  if  ■  Pi"  ii]  i,| 
if  ij  ^  ^ 


a  = 


il 

1  20i+'2)  iz  ^  1  ij 

2 

[  3  /,  12  2j,  12. 

(2.6-47a) 


Ij  ^  Pi"  ij  -  Pj’  ii  -  a  it  h  ih  +  h) 

a  i}  il 


{2.6-Alb) 


^  =  i\  Pi'-iiPi"  +  ah  il 

■2  1 
a  if  li 


(2.6-47C) 
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Determination  of  smoothed  fluxes  in  region  1 .  The  fluxes  that  we  actually  record  for 
the  grid  points  within  region  1  should  be  related  to  the  simple  cubic  formula  (2.6-44<3)  (and 
(2.6-44^7))  by  a  smoothing  integration.  We  interpret  P’'(IJ)  as  given  in  (2.6-44a)  and 
(2.6-446)  as  indicating  the  local  precipitation  rate  at  grid  point  (7,7).  Define  P/XU)  to  be 
the  integral  of  the  continuous  function  P'XiJ)  from  /  -  0.5  to  7  +  0.5.  Assuming  that 
IMIN  -  it  >  0.5,  which  should  normally  be  the  case,  then,  for  grid  points  in  the  range 
i/+0.5  <  7  <  IMIN,  we  have 


where  /+  =  7  -t-  0.5  -  it  and  i.  =  7  -  0.5  - 1; .  For  the  7-value  that  is  centered  on  the 
segment  of  the  /-grid  line  that  is  cut  by  ellipse  #1,  i.e.,  the  7-value  that  lies  closest  to  ellipse 
#1,  we  have  to  cut  off  the  integral  at  i/  rather  than  7  -  0.5,  and  we  get 


P  " 
‘  s 


a 


il  + 


(2.6-49) 


In  the  unlikely  event  that  IMIN  -  i/  <  0.5,  then  we  should  use  formula  (2.6-49)  for  7  = 
IMIN,  and  /*,"  =  0  for  smaller  7-values.  In  any  case,  I  propose  that  Ps'  -  0  for  7  <  it 
-0.5. 


Emergency  procedure  in  the  case  where  7^"  <  0  for  some  grid  point.  In  this  case, 
which  should  occur  extremely  rarely,  we  simply  make  Ps"  the  same  for  all  7  in  the  range  // 
-  0.5  <  7  <  IMIN,  and  make  the  total  satisfy 

X  PsVrT)  +  05  PsVMINJ)  =  7>r(7) 

ir0.5<I<IM[N 
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2.7.  Limitations  of  the  Magnetospheric  Specification  Model. 

The  MSM  is  the  first  effort  at  a  large-scale  computational  model  of  the  Earth's 
magnetosphere  for  operational  use.  As  indicated  by  Ae  discussion  of  test  results  in 
Sections  3  and  4  of  this  report,  this  effort  has  been  a  success.  With  regard  to  its  primary 
goal  of  specifying  fluxes  of  <  100  keV  electrons  in  the  geostationary-orbit  region,  the 
MSM  does  a  very  good  job.  'The  MSM  calculates  these  geostationary-orbit-region 
electrons  within  a  comprehensive  framework  that  also  provides  flux  values  throughout  the 
middle  magnetosphere,  as  well  as  magnetospheric  inputs  to  the  ionosphere  and 
thermosphere. 

At  the  same  time,  the  MSM  has  not  reached  the  ultimate  goal  of  being  capable  of 
specifying  all  important  physical  parameters  at  all  locations  in  the  magnetosphere  with  high 
confidence  of  high  accuracy  for  all  conceivable  conditions.  The  limitations  result  from 
three  basic  sources:  (i)  Our  knowledge  of  the  Earth’s  magnetosphere  is  far  from  complete; 
(ii)  'The  input  data  that  is  expected  to  be  available  to  the  Space  Forecast  Center  is  far  from 
ideal;  (iii)  Only  a  very  limited  amount  of  observational  data  has  been  available  to  test  the 
model. 

With  regard  to  the  MSM's  primary  goal  of  correctly  specifying  fluxes  of  1-100  keV 
electrons  in  the  geostationary-orbit  region,  the  essential  limitations  are: 

A 1 .  The  MSM  does  not  correctly  represent  flux  dropouts  in  the  >  40  keV  electrons,  and  it 
sometimes  specifies  high  fluxes  when  only  average  levels  are  observed.  The  operational 
impact  of  these  deficiencies  is  reduced  by  the  fact  that  the  MSM  automatically  monitors  and 
reports  on  its  consistency  with  real-time  observations  from  geostationary  spacecraft. 

A2.  Testing  of  the  model  has  been  limited  by  the  available  observational  database. 
Although  we  expect  the  model  to  perform  well  for  geostationary-orbit  electrons  below  35 
keV,  no  spacecr^t  data  have  been  available  to  test  the  MSM  in  this  energy  range. 

A3.  The  model  has  been  tested  for  quiet  periods  and  substorms,  up  through  a  moderately 
strong  magnetic  storm.  It  has  not  been  possible  yet  to  assemble  a  test  data  stream  for  a 
very  large  magnetic  storm. 

A4.  Although  the  MSM  has  been  tested  for  two  events  that  occurred  nine  years  apart,  it 
has  not  been  tested  for  all  phases  of  the  solar  cycle. 

Two  limitations  that  follow  from  the  general  formulation  of  the  MSM  are  the 
following: 

B 1 .  Its  capabilities  for  forecasting  are  very  limited  unless  an  upstream  solar-wind  monitor 
is  available. 

B2.  It  does  not  keep  track  of  particle  pitch-angle  distributions. 

With  regard  to  the  broader  goal  of  specifying  all  important  magnetospheric  parameters 
throughout  the  middle  magnetosphere,  the  capabilities  of  the  MSM  have  generally  not  yet 
been  as  well  tested  as  we  would  like.  Nevertheless,  we  believe  that  the  MSM  makes  useful 
estimates  of  electric  fields,  magnetic  fields,  and  particle  fluxes  (up  to  ~  3(X)  keV)  for  the 
global  ionosphere  and  for  the  region  3<L<10  of  the  magnetosphere.  The  following 
specific  limitations  apply: 

(Jl.  No  flux  estimates  are  currently  made  for  electrons  above  about  300  keV  at 
geostationary  orbit  (roughly  (3(X)  keV)(L/6.6)'*/3  for  other  L-shells).  H.  C.  Koons  and  D. 
J.  Gomey  have  recently  provided  us  with  their  algorithm  for  specifying  and  predicting 
>  3  MeV  electrons  for  the  geostationary-orbit  region.  However,  that  code  was  delivered 
after  the  April  1  end  of  the  present  contract,  and  we  consequently  have  not  yet  been  able  to 
install  and  test  the  Koons-Gomey  algorithm  within  the  framework  of  the  MSM.  A  new 
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version  of  the  MSM  that  includes  the  Koons-Gomey  algorithm  will  be  delivered  soon  after 
July  1,  1990. 

C2.  For  100-300  keV  electrons  and  50-300  keV  ions,  the  MSM  calculates  fluxes  by 
scaling  an  internally  computed  set  of  base-level  fluxes  to  fit  the  real-time-observed 
geostationary  fluxes. 

C3.  The  accuracy  of  the  MSM  for  electrons  <  100  keV  (6.6/L)*/3  becomes  increasingly 
uncertain  as  L  decreases  far  below  6.6,  because  our  test  data  sets  contained  no  spacecraft 
observations  inside  L  =  6.6. 

The  present  uncertainties  concerning  model  performance  should  decrease  rapidly  over 
the  next  several  years.  A  data  stream  is  being  collected  for  the  great  ma^etic  storm  of 
March  1989,  and  MSM  tests  on  that  event  will  begin  soon:  this  should  eliminate  limitation 
A3.  For  the  period  beginning  in  the  Fall  of  1989,  additional  particle  detectors  are  being 
placed  on  DoD  geostationary-orbit  spacecraft:  these  new  detectors,  developed  at  Los 
Alamos,  monitor  fluxes  of  low-energy  electrons  and  ions;  we  expect  to  assemble  data  sets 
to  test  the  MSM,  including  those  data,  within  the  next  year,  thus  elimirating  limitation  A2. 
The  CRRES-SPACERAD  spacecraft,  scheduled  for  launch  in  July  1990,  will  provide  a 
rich  source  of  data  for  L  <  6.6;  thus  within  the  next  1-2  years,  limitation  C3  will  be 
eliminated.  With  this  continued  testing,  limitation  A4  will  gradually  be  eliminated.  The 
flexible,  modular  structure  of  the  MSM  will  allow  us  to  carty  out  simple  "fixes"  by 
changing  very  few  lines  of  code,  if  the  extended  testing  indicates  the  need  for  such 
changes. 
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2.8  MSM  Output 

The  primaiy  output  of  the  MSM  is  the  parti''le  flux,  for  each  energy,  and  species,  on 
each  giid  point.  The  fluxes  are  labeled  by  the  indices  I,  J,  and  K  where  I  is  the  index 
specifying  the  geomagnetic  colatitude  and  J  specifies  the  magnetic  local  time  of  the  grid 
point.  K  is  the  invariant  energy  index. 

Fluxes  are  also  specified  on  the  model  boundary. 

The  X.Y.  and  Z.  position  of  each  grid  point  on  the  equatorial  reference  surface  is 
provided. 

For  the  precipitating  particles,  the  energy  flux  and  average  energy  of  the  electrons  and 
ions  is  given  at  each  grid  point. 

A  time  tag  is  provided  that  gives  the  time  of  each  output  record. 

Numerous  other  parameters  are  included  in  the  output  record.  These  include  the  input 
data  parameters  actually  used  by  the  model,  as  well  as  a  large  number  of  internal  model 
parameters  that  can  be  used  to  restart  the  model  or  for  diagnostic  purposes. 

Table  2.8  gives  the  entire  list  of  output  parameters  together  with  their  program  names 
and  units. 

The  parameter  ERSHFT  specifies  the  average  difference  between  the  model's 
prediction  for  the  log  of  the  geostationary  flux  and  the  average  observed  value  for  the  last 
15  minutes,  for  each  MSM  energy  channel.  This  could,  for  example,  provide  a  warning  of 
an  unexpected  situation  where  the  MSM  is  substantially  underestimating  the  flux. 
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Table  2.8 


MSM  Output  Parameters 


PARAMETER 


NAME 


UNITS 


GRID  POINT 

VM(IJ) 

1 

Re/nT 

PART.  FLUX  <S)  GRID  POINTS  U  AND  ENERGY  K 


PART.  FLUX  ON  BOUND.  J  &  ENERGY  K 


X  POSITION  OF  GRID  POINT  U  ON  EQUAT.  REF.  SURF. 


Y  POSITION  OF  GRID  POINT  IJ  ON  EQUAT.  REF.  SURF. 


Z  POSITION  OF  GRID  POINT  I,J  ON  EQUAT.  REF.  SURF. 


ENERGY  FLUX  OF  ELECTRONS  @  IONOSPHERE  G.P.  I.J 


AVERAGE  ENERGY  OF  ELECTRONS  @  IONOSPHERE  G.P.  I J 


ENERGY  FLUX  OF  IONS 


AVERAGE  ENERGY  OF  IONS  @  IONOSPHERE  G.P.IJ 


POTENTIAL  @  IONOSPHERIC  G.P.  I J 


MAGNETIC  COLATITUDE  OF  IONOSPHERIC  G.P.  U 


MAGNETIC  LOCAL  TIME  OF  IONOSPHERIC  G.P.  U 


FLUX  (IJ.K 


FLXBND(JJ< 


MAGNETOPAUSE  STANDOFF  DISTANCE 


DIPOLE  TILT  ANGLE 


TAIL  CURRENT  SHEET  DENSITY 


RATIO  TAEL  I O  INFINITY  TO  TAIL  I  &  INNER  EDGE 


Y  VARIATION  OF  CROSS-TAIL  CURRENT 


CURRENT  SHEET  THICKNESS  PARAMETER 


DEGREE  OF  COLLAPSE  MIDNIGHT  REGION  TAIL 


RANGE  IN  X  AFFECTED  BY  TAIL  COLLAPSE 


PARM.  DETER.  THE  Y-DEP  OF  COLLAPSE  CURRENT 


-B  PERT.  (®  EARTH  CENTER  DUE  TO  WEST  RING  I 


+B  PERT.  @  EARTH  CENTER  DUE  TO  EAST  RING  I 


RADIUS  OF  THE  WESTWARD  RING  CURRENT 


RADIUS  OF  THE  EASTWARD  RING  CURRENT 


MAGNETIC  FIELD  STRENGTH  ON  EQUAT  REF  SURFACE 


POLAR  CAP  POTENTIAL  DROP 


POLAR  CAP  PATTERN  TYPE 


RATE  OF  MOTION  OF  THE  LOW  LAT.  AURORAL  BOUNDARY 


BOUNDARY  2  &  3  ELLIPSE  SPECmCATION  PARAMETER 


BOUNDARY  2  &  3  ELLIPSE  SPECIFICATION  PARAMETER 


BOUNDARY  2  &  3  ELLIPSE  SPECIFICATION  PARAMETER 


BOUNDARY  2  &  3  ELLIPSE  SPECIHCATION  PARAMETER 


ELECTRIC  POTENTIAL  N.  IONOSPHERE 


ELECTRIC  POTENTIAL  S.  IONOSPHERE 


BOUNDARY  LOCATION  IN  OUR  GRID 


LOGARITHMIC  ERROR  SHIFT 


PROCESSED  INPUT  ARRAY 


TIME  TAG  OF  RECORD  TIME 


EAVGaJ.2) 


vaj 


COLATCU 


ALOCT(U) 


STFD 


TILT 


HJNEAR 


HJFRAC 


DDY 


DD 


DELXCL 


DYC 


BRN 


BRP 


RHRN 


RHRP 


PCP 


IPATT 


DXfL) 


DY(L) 


VNORTHfU) 


VSOUTHfU 


BNDLOCU) 


ERSHFUK) 


AUGPAR 


TIME 


KEV 


VOLTS 


RADIANS 


RADIANS  east 
from  noon 


Re 


DEGREES 


N/A 


N/A 


Re 


Re 


N/A 


Re 


Re 


NANOTESLA 


NANOTESLA 


Re 


Re 


NANOTESLA 


KV 


N/A 


DEGREES/HR 


DEGREES 


DEGREES 


DEGREES 


DEGREES 


VOLTS 


VOLTS 


I  GRID  UNITS 


VARIOUS  UNITS 


YR:DAY  SEC 


MSM  SIMPLIFIED  FLOW  CHART 


Figure  2.1 


sfc-6/21/90 


BlHil 


Figure  2.2 


Data  from 
Previous 
MStM  Runs^ 


Real-Time 

Satellite 


Real-Time 

Ground-Based 

Data 


6/21/90 


Environmental  Datalbase 
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Latitude  of  Equatorward  Edge 


Figure  2.3-1.  Formulas  relating  the  equatorial  distance  of  the  inner  edge  of  the  plasma 
sheet  at  local  midnight  to  the  latitude  of  the  equatorward  edge  of  the  aurora  at  local 
midnight.  The  curves  marked  r-high  and  r-low  represent  and  upper  lower  envelope  of  a  set 
of  estimates  derived  from  the  papers  by  [Horwitz,  1986  #9;  Kivelson,  1976  #8; 
Gussenhoven,  1983  #3].  The  curve  marked  raipoie  represents  mapping  in  a  simple  dipole 
field,  and  the  curve  marked  r-our  formula  represents  the  result  of  applying  equation 
(2.3.1). 
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Polar  cap 


Figure  2.4- 1 .  Division  of  the  ionosphere  into  four  regions  for  the  electric-field  model. 
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Figure  2.4-1,  The  assumed  variation  of  with  <])  in  the  low-latitude  region  of  the  electric- 
field  model,  with  the  minimum  value  normalized  to  -1,  and  the  normalized  indefinite 
integral  of  ,  as  indicated  in  equation  (2.4-20). 
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Figure  2.4-2.  Plots  of  meridional  electric  fields  based  on  equations  (2.4-29)-(2.4-31). 
The  top  plot  pertains  to  A0  =  0.5,  the  lower  plot  to  A0  =  2.0. 
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Figure  2.4-4.  Equipotential  diagram  computed  for  0045  UT  on  April  22,  1988.  The  view 
is  of  the  northern  ionosphere,  with  noon  toward  the  top  of  the  diagram.  There  is  a  6  kV 
potential  difference  between  adjacent  equipotentials. 
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Figure  2.4-5.  Equipotential  diagram  computed  for  0300  UT  on  April  22, 1988.  The  view 
is  of  the  northern  ionosphere,  with  noon  toward  the  top  of  the  diagram.  There  is  a  6  kV 
potential  difference  between  adjacent  equipotentials. 
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Figure  2.5-1.  Observational-average  and  kappa-function  geosynchronous-orbit  electron 
fluxes  for  quiet,  pre-substorm  conditions.  The  sources  of  observational  data  are  given  in 
Table  2.5-3. 
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Figure  2.5-2.  Flux  values  for  L  =  3  and  4  derived  from  NASA  Model  AE7,  and  our  fitting 
kappa  functions. 
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Figure  2.5-3.  Fluxes  Jkp  of  plasma-sheet  hydrogen  and  oxygen  ions. 
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Figure  2.5-4.  Hydrogen  and  oxygen  fluxes,  as  represented  in  jup,  for  geosyncnronous 
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H+  Fluxes  fo#Cp  =  5+  Storm  Time 


0+  Fluxes  foiKp  =  5+  Storm  Time 


Figure  2.5-6.  Fit  of  model  fluxes  (m)  to  observational  data  (d)  for  an  inbound  AMPTE 
pass  that  occurred  early  in  the  day  September  5,  1984.  The  data  source  was  Figure  1 2a  of 
[Gloeckler,  1987  #23).  Kp  at  this  time  was  5+. 


14 


H+  Fluxes  forKp  =  9-  Storm  Time 
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Figure  2.5-7.  Fit  of  model  fluxes  (m)  to  observational  data  (d)  for  an  inbound  AMPTE 
pass  that  occurred  early  on  February  9,  1986.  The  data  source  was  Figure  8a  of 
[Hamilton,  1988  #24].  fCp  at  this  time  was  9-. 
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Figure  2.6-3.  Quarter-bounce  times  for  H'*’  and  0+  as  functions  of  energy  for  the  indicated 
L-values  (r^^).  Solid  curves  denote  times  for  pitch  angles  very  near  the  loss  cone,  while 
dashed  curves  are  for  particles  mirroring  very  near  the  equatorial  plane.  The  loss  cone  has 
been  taken  to  be  defined  by  the  planet  surface  (dipole  magnetic  field). 
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Figure  2.6-4.  O-eq)  results  for  r^a  =  3.0  Re  and  6.6323  Re  as  a  function  of 

equatorial  pitch  angle.  The  loss  cone  is  defined  by  the  planet  surface.  The  exobase 
temperature  in  these  calculations  has  been  taken  to  be  1000°  K  and  the  exobase  density  8.0 
X  10^  cm-T 
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Figure  2.6-5.  Stationary-case  decay  (erosion)  of  ring  current  ions  is  illustrated  as  a 
function  of  particle  energy  by  plotting  the  pitch-angle  inte^ated  KDF  -  denoted  n{E,t)  -- 
at  six-hour  intervals,  normalized  by  n(£,0).  The  loss  cone  is  defined  by  the  planet  surface, 
the  exobase  temperature  is  1(X)0°  K,  and  the  exobase  density  is  8.0  x  10^  cm-3.  (a)  H+  at 
3.0  Re. 
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Figure  2.6-5b.  O'*’  at  3.0  Re 
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Figure  2.6- 5c.  at  6.6323  Re-  (Note  the  change  in  vertical  scale.) 
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Figure  2.6- 5d.  0+  at  6.6323  Re- 
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3.  TEST  EVENTS:  MARCH  22,  1979  AND  APRIL  21-23,  1988 

3.0.  Introduction 

The  Magnetospheric  Specification  Model  has  been  tested  against  two  major 
geomagnetic  storm  events  separated  in  time  by  nine  years.  The  first  test  event  occured  on 
March  22,  1979,  about  1  year  before  the  peak  of  solar  cycle  21.  This  event  was  chosen 
because  it  is  a  well-documented  event  (the  CDAW  6  event)  and  a  large  body  of  data  from  a 
variety  of  sources  was  readily  available  to  compare  with  the  model  output,  however, 
DMSP  ion  drift  meter  data  were  not  available  for  this  event.  The  second  event  occured 
April  21-23, 1988,  about  2  years  before  the  peak  of  solar  cycle  22.  The  second  event  was 
chosen  because  GL  had  already  begun  an  investigation  of  the  event  and  DMSP  data  were 
available.  In  both  cases,  the  choice  of  events  was  dictated  by  the  availability  of  test  data. 
Both  events  represent  important  geomagnetic  storms. 

3.1.  The  March  22,  1979  Event 

The  March  22,  1979  event  was  used  about  mid-way  through  the  contract.  This  was 
before  the  MSM  was  developed  to  the  point  of  fully  automated  operation.  Model  time  steps 
and  input  data  were  hand  tailored  for  the  event.  In  this  respect,  it  was  not  a  test  of  the 
MSM  in  it's  delivered  version.  Moreover,  the  lack  of  DMSP  ion  drift  meter  data  limited  the 
electric  field  model.  For  these  reasons  we  will  not  discuss  this  event  further  except  to  note 
that  the  success  of  the  model  with  this  event  provided  confidence  that  the  basic  design  of 
the  MSM  was  valid  and  formed  the  basis  of  the  Critical  Design  Review. 

The  detailed  results  of  the  MSM  using  this  event  can  be  found  in  Quanerly  Status 
Report  nos.  6  and  7  and  in  Bonnie  Hausman's  M.S.  thesis.  This  thesis,  which  contains  a 
good  description  of  the  preliminary  version  of  the  MSM,  is  included  in  the  Appendix  of 
this  repon. 

3.2.  The  April  21-23,  1988  Event 

3.2.1  Geophysical  Conditions  for  the  Event 

As  can  be  seen  from  figure  3.1,  which  shows  five  of  the  geophysical  parameters 
available  as  input  to  the  MSM,  the  April  21-23, 1988  consisted  of  a  large  initial  substorm,  a 
period  of  strong  but  unsteady  convection,  a  quiet  period,  and  then  a  second  large  storm 
almost  exactly  one  day  after  ^e  start  of  the  first  substorm.  The  standoff  distance  shown  in 
figure  3.1  was  derived  from  solar  wind  parameters  shown  in  figure  4.3  and  4.4.  The 
frame  labeled  Equatorward  Edge  is  the  low-latitude  boundary  of  the  auroral  zone  at  the 
midnight  meridian  as  derived  from  the  electron  detectors  on  the  DMSP  satellites.  These 
data  were  provided  by  Bill  Denig  of  GL.  Figure  3.2  shows  the  fluxes  obtained  by  three 
geostationary  spacecr^t  during  the  event.  These  data  are  from  detectors  designed  and  built 
by  the  Los  Alamos  National  Laboratory.  We  had  help  interpreting  the  detector  sensitivities 
and  geometric  factors  from  Tom  Cayton  and  the  staff  at  LANL. 

From  figures  3.1  and  3.2,  we  see  that  the  onset  of  the  first  substorm  begins  shortly 
before  the  start  of  day  113  (April  22nd).  It  rises  to  a  peak  by  0200  hrs  on  1 13  and  remains 
intense  until  about  12(X)  hrs.  Dst  reached  it's  most  negative  value  of  -120  nanotesla  around 
1000  hrs  on  day  113,  indicating  that  the  storm-time  ring  current  peaked  at  that  time.  The 
second  large  substorm  occurs  part  way  into  the  recovery  phase  of  this  storm,  at  the  start  of 
day  114. 

3.2.2.  Comparison  With  the  MSM  Output 


Section  3,  Page  2 


In  order  to  study  the  MSM  predicted  fluxes  in  the  equatorial  plane  and  to  compare 
these  fluxes  with  the  spacecraft  data  we  have  prepared  an  animat^  video  of  the  MSM 
output  for  this  event.  In  this  video,  the  MSM  fluxes  are  color  coded  and  the  model  output 
after  each  15  minute  time  step  is  speeded  up  under  animation.  The  location  of  each  satellite 
is  superimposed  on  the  magnetospheric  equatorial  plane  in  the  model  output  display  so  it  is 
easy  to  see  where  each  satellite  is  throughout  the  storm. 

In  the  video,  as  time  progresses,  the  satellites  advance  around  their  orbits,  and  the 
storm  front  (the  inner  edge  of  the  plasma  sheet)  is  seen  to  intensify,  move  toward  the  Earth, 
and  cross  the  orbit  of  the  satellites.  The  format  of  a  typical  video  frame  is  shown  in  figure 
3.3.  The  MSM  output  is  shown  in  the  upper  left-hand  side  of  the  frame.  In  the  model 
frame,  the  sun  is  to  the  left  and  the  dawn  meridian  is  at  the  top.  The  fluxes  plotted  are  40 
KeV  electrons.  The  30  -44  KeV  electron  fluxes  reported  by  satellite  1  are  shown  beside  the 
MSM  output  on  the  upper  right,  and  four  of  the  relevant  model  input  parameters  are  shown 
at  the  bottom  of  the  frame.  In  both  the  satellite  flux  frame  and  the  parameter  frame  below, 
the  time  being  portrayed  by  the  model  output  is  shown  as  a  dashed  vertical  line. 

The  video  illustrates  very  dramatically  the  success  of  the  MSM  at  predicting  the  fluxes 
observed  at  satellite  #1.  Unfortunately  it  is  not  practical  to  include  the  video  with  this 
report  so  we  will  instead  walk  through  a  sequence  of  seven  frames,  consisting  of  figures 
3.3  through  3.9,  which  illustrate  interesting  times  in  the  event. 

Figure  3.3  illustrates  a  point  in  the  growth  phase  of  the  first  substorm.  Notice  that 
satellite  1  is  midway  between  the  dusk  and  midnight  meridians  and,  according  to  the 
model,  and  the  observations,  has  not  yet  seen  a  flux  enhancement.  The  storm  front,  as 
indicated  by  the  orange  region,  has  not  yet  advanced  to  the  geostationary  orbit.  The  model 
and  the  observations  agree  at  this  point. 

Figure  3.4  shows  the  situation  a  short  time  later,  at  the  substorm  onset  at  the  satellite. 
Note  that  the  MSM  now  shows  satellite  1  right  at  the  edge  of  the  storm  front  (large  flux 
gradient)  as  indicated  by  the  proximity  of  the  orange  and  deep  orange  regions  and  the 
satellite  flux  graph  shows  a  large  flux  enhancement  about  to  begin. 

Figure  3.5  now  shows  satellite  1  fully  engulfed  in  the  storm-enhanced  plasma  as  seen 
in  both  the  MSM  output  and  the  satellite  observed  fluxes. 

Figure  3.6  shows  the  peak  of  the  storm  as  seen  by  satellite  1  and  as  predicted 
perfectly  for  that  satellite  by  the  model. 

In  figure  3.7,  the  satellite  has  moved  to  the  dawn  meridian  and  is  heading  toward  the 
region  of  decreasing  flux  values  on  the  day  side  of  ihe  magnetosphere.  The  satellite 
observations  already  show  a  low  flux  value  and  so  the  model  is  somewhat  behind. 

In  figure  3.8,  nearly  one  day  has  elapsed  since  the  start  of  the  sequence,  and  we  find 
satellite  1  about  to  enter  a  storm  front  again,  this  time  for  the  second  substorm.  Again,  the 
model  agrees  well  with  the  observations. 

Finally,  in  figure  3.9,  we  are  in  the  recovery  phase  of  the  second  substorm.  The 
storm  front  has  receded  outside  the  geostationary  orbit.  Low  fluxes  are  predicted  by  the 
model  and  observed  by  satellite  1. 
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We  believe  that  this  sequence  of  figures  provides  proof  that  the  MSM  does  a  good  job 
in  modeling  the  timing  of  flux  enhancements.  In  section  4,  we  discuss  the  accuracy  of  the 
MSM  from  a  statistical  standpoint. 

Comparison  of  the  MSM  with  the  other  two  satellites,  for  this  event,  is  possible  by 
examination  of  figures  4.9  through  4.14,  in  section  4.  These  figures  show  the  detailed 
time  comparison  between  the  MSM  output  fluxes,  interpolated  to  the  three  satellite  orbit 
positions,  for  40  KeV  and  65  KeV  electrons.  These  figures  also  have  the  Garrett  model 
output  displayed.  Section  4  contains  a  discussion  of  these  figures  and  conclusions  on  the 
overall  accuracy  of  the  MSM. 
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4.  VALIDATION  STATISTICS 

In  order  to  obtain  statistics  to  determine  the  accuracy  of  the  MSM,  the  model  output 
was  compared  with  electron  flux  data  from  three  geostationary  orbit  satellites  supplied  by 
the  Air  Force.  No  geostationary  satellite  flux  data  were  used  as  inputs  to  the  model.  The 
three  satellites  were  well  spaced  in  local  time  and  therefore  represent  a  good  test  of  the 
model  accuracy  at  all  local  times.  The  event  used  was  the  large  magnetic  storm  of  April  21- 
23,  1988  discussed  in  section  3.2.  The  following  sections  describe  the  results  of  this 
study. 


4. 1 .  Model  Inputs 

The  input  data  to  the  MSM  used  for  the  test  event  are  Kp,  Dst,  solar  wind  velocity, 
solar  wind  density,  the  polar  cap  potential  drop,  the  polar  cap  potential  pattern  type  and  the 
equatorward  edge  of  the  auroral  zone  (see  Table  4.1).  The  solar  wind  velocity  and  density 
were  used  to  calculate  the  standoff  distance,  which  is  the  distance  from  the  center  of  the 
earth  to  the  magnetopause  at  the  subsolar-point.  The  standoff  distance,  Dst  and  the 
equatorward  edge  of  the  auroral  zone  are  used  by  the  B-field  model  in  the  MSM.  The  E- 
field  model  uses  the  polar  cap  potential,  the  polar  cap  potential  pattern,  the  equatorw'ard 
edge  of  diffuse  aurora  and  its  time  derivative.  Initial-  and  boundary-condition  fluxes  are 
based  on  Kp. 


Table  4.1. 


The  MSM  Input  Parameters 


Input 

Base 

Figure 

Kp 

Ground  observation 

Figure  4.1 

Dst 

Ground  observation 

Figure  4.2 

Solar  wind 
velocity 

Satellite  observation 
NSSDC 

Figure  4.3 

Solar  wind 
density 

Satellite  observation 
NSSDC 

Figure  4.4 

Polar  cap 

potential 

drop 

Satellite  observation 
DMSP  data 

GL/UTD 

Figure  4.5 

Polar  cap 

potential 

pattern 

Satellite  observation 
DMSP  data 

GL/UTD 

Figure  4.6 

Equatorward 
edge  of  the 
auroral  zone 

Satellite  observation 
DMSP  data 

GL 

Figure  4.7 

NSSDC;  The  National  Space  Science  Data  Center 
GL;  Geophysical  Laboratory 
UTD;  University  of  Texas  at  Dallas 
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4.2.  Model  Outputs 
4.2.1  The  Interpolation 

In  order  to  compare  the  MSM  output  with  the  satellite  ob'^ervation,  the  MSM 
predicted  fluxes  had  to  be  determined  at  the  three  Air-Force  satellite  local  times  and 
locations.  Since  the  output  time  intervals  of  the  MSM  nd  the  satellite  data  times  are 
generally  different,  and  since  the  MSM  grid  poi..  ^  are  not  necessarily  at  the 
geosynchronous  orbit  local  times  of  the  satellites,  the  model  output  fluxes  had  to  be 
interpolated  in  time  and  space.  This  interpolation  was  accomplished  as  follows: 

(1)  First,  the  spatial  interpxalation  was  performed  both  for  ti  and  t2  where 


ti  <  tsatellite  ^2 

and  ti,  t2  are  the  times  of  the  closest  MSM  output.  This  spatial  interpolation  was  linearly 
calculated  in  terms  of  the  fluxes  at  the  model  grid  points  as  shown  in  figure  4.8. 


(2)  Then  the  two  model  fluxes  J(ti)  and  J(t2 )  were  linearly  interpolated  in  time: 
J(at  the  satellite)  =  J(ti)*(t2  -tsatellite)/(t2  -tl)  +  J(t2  )*(tsatemte-tl)/(t2  -tl). 


In  summary,  this  is  a  3-dimensional  (2D-spatial  and  ID-temporal)  linear 
interpolation.  In  this  evaluation,  the  magnetic  and  the  geographic  equatorial  planes  are 
assumed  to  be  the  same.  And  the  quantity  (flux  volume)''(-2/3)  at  the  geosynchronous 
orbit  is  taken  to  be  approximately  7.0. 


4.2.2  Results 

The  electron  fluxes  at  geosynchronous  orbit  were  computed  for  two  differential 
energy  channels,  30-44keV  and  44-64keV,  using  the  data  supplied  by  the  Air  Force.  The 
differential  channel  passbands  of  the  detectors  are  much  larger  then  these  widths  and  so  the 
fluxes  in  adjacent  channels  were  subtracted  to  form  the  smaller  differential  channels.  In 
consultation  with  Tom  Cayton  of  LANL,  a  standard  correction  was  added  for  detector 
dead-time.  At  the  same  time,  the  low  st  channel  flux  was  increased  by  a  factor  of  2  and  the 
next  lowest  by  1.4  to  achieve  agreement  with  ATS  data  [Garrett,  Private  Communication]. 
These  corrections  are  discussed  in  detail  in  Quarterly  Report  no.  1 1 . 

The  MSM  output  fluxes  and  the  corresponding  satellite  fluxes  are  plotted  in  figures 
4.9  to  4.14.  In  these  figures,  the  thick  solid  lines  denotes  the  MSM  output,  the  solid  lines 
with  square  markings  are  the  satellite  measurements.  Also  shown  as  the  solid  thin  lines  are 
the  output  from  a  model  of  the  average  geostationary  electron  fluxes  prepared  by  Garrett . 
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4.3.  The  Garrett  model 

Through  the  courtesy  of  H.B. Garrett,  Geophysics  Laboratory  (currently  at  the  Jet 
Propulsion  Laboratory),  we  obtained  a  statistical  model  that  computes  average  electron 
fluxes  at  arbitrary  loctil  time  for  Kp  values  up  to  6  at  the  geostationary  orbit.  This  model  is 
described  in  detail  in  section  2.2. 2.4.  The  Garrett  model  was  run  for  the  April  21-23,  1988 
event.  The  model  generates  geostationary  electron  fluxes  with  the  input  of  the  satellite  local 
time  and  Kp. 


4.4.  Evaluation  of  Accuracy 

For  the  e%  aluation  of  the  accuracy  of  a  model,  it  would  seem  to  be  desirable  to  have  a 
single  parameter,  or  accuracy  index,  which  can  be  taken  as  a  indicator  of  the  extent  to 
which  the  model  follows  the  observations.  In  fact,  as  we  shall  see,  such  an  approach  is 
risky  when  the  goal  is  to  build  a  model  which  accurately  predicts  occurrences  of  worst-case 
or  adverse  storm  conditions.  Nonetheless,  we  have  chosen  for  such  an  index  the  root- 
mean-square  of  the  log  of  the  ratio  of  the  model  to  the  satellite  fluxes  averaged  over  the 
most  active  two  days  of  the  test  event.  It  is  essential  to  use  the  logarithm  because  of  the 
large  dynamic  range  of  the  fluxes,  nearly  four  orders  of  magnitude. 


log  10 


Jmodel 

Jsatellite 


This  parameter  was  calculated  for  each  of  the  two  differentia!  energy  bands  mentioned 
above  and  for  all  tl.  'ee  satellites.  The  results  are  shown  in  Tables  4.2  and  4.3. 


Table  4.2.  The  rms  of  errors  (40keV) 


model 

Satellite- 1 

Satellite-2 

Satellite-3 

MSM 

0.7350 

0.5949 

0.7434 

Garrett 

1.2062 

1.1683 

1.2817 

Table  4.3.  TTie  rms  of  errors  (65keV) 

model 

Satellite- 1 

SateIlite-2 

SateIlite-3 

MSM 

0.9238 

0.7630 

0.6380 

0.5328 

0.4758 

It  is  seen  that  the  best  agreement  is  with  satellite  -2  at  the  higher  energy.  For  this  case  the 
difference  would  be  a  factor  of  3.3.  The  worst  case  is  found  for  satellite  -1  for  the  same 
energy  range,  a  factor  of  8.4. 
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4.5  Conclusions 

Based  on  the  time-dependent  data  presented  in  this  section  and  section  3.2, 
throughout  the  run,  the  Magnetospheric  Specification  Model  was  able  to  predict  storm 
related  electron  enhancements  in  good  dme  agreement  with  geosynchronous  satellite  data. 

Based  on  the  accuracy  index  reported  in  section  4.4,  the  MSM  prediction  was  very 
much  better  than  the  statistical  Garrett  model  for  40keV  electron  fluxes.  This  is  panly 
because  the  Garrett  average  fluxes  were  too  high  for  this  energy  range.  However,  for  the 
higher  energy  channel,  the  index  showed  that  the  MSM  predictions  were  worse  than  the 
Garrett  model.  At  the  same  time,  the  MSM  fluxes  lollow^  enhancements  in  the  fluxes  far 
better  then  the  Garrett  model.  T^e  Garrett  model  fails  completely  in  predicting  the  high 
fluxes.  This  illustrates  the  inappropriateness  of  a  single  parameter  index  for  the  purpose  of 
assessing  model  accuracy. 

The  MSM  never  failed  to  predict  high  electron  fluxes  when  they  were  observed. 
There  were,  however,  some  instances  when  the  satellites  reported  sudden  decreases,  or 
flux  dropouts,  not  predicted  by  the  MSM.  We  tentatively  attribute  these  instances  to  two 
possible  situations.  In  some  cases,  the  plasma  sheet  may  have  been  very  thin,  as  reported 
to  be  possible  by  Sergeev  et  al.,  [19901,  so  that  the  spacecraft  may  have  been  beyond  and 
above  or  below  the  plasma  sheet  and  hence  on  open  field  lines  or  on  field  lines  that  extend 
very  far  out  into  the  plasma  sheet  and  therefore  contain  lower  plasma  fluxes.  Our 
magnetic-field  model  is  not  streched  enough  to  represent  these  very  thin  plasma  sheet 
configurations.  In  other  cases,  the  electric-field  model  may  have  resulted  in  the  particles 
being  injected  too  deeply  into  the  magnetosphere.  The  sunward-flow  region  of  the  electric- 
field  model  tends  to  be  more  confined  in  latitude  than  is  consistent  with  observations, 
which  would  account  for  this  extra-deep  injection.  There  has  not  been  sufficient  time  in  the 
contract  to  correct  these  late-surfacing  and  rather  minor  deficiencies,  however,  a  B-field 
model  with  a  thinner  current  sheet  is  being  tested. 

Near  0200  UT  of  day  1 14,  satellite-2  observed  an  average  quiet-time  electron  flux 
level,  but  the  MSM  predicted  a  flux  dropout.  At  the  same  time  that  the  dropout  was 
predicted  by  the  MSM,  a  strong  enhancement  of  the  polar  cap  potential  was  observed  by 
the  DMSP  satellite  (see  figure  4.5).  The  MSM  dropout  represented  the  transport  of  low- 
density  flux  tubes,  which  had  previously  been  on  trapped  orbits  at  the  geostationary  orbit, 
to  untrapped  orbits  which  intersect  the  magnetopause.  The  fact  that  Spacecraft-2  did  not 
experience  any  dropout  may  indicate  that,  either  the  model  overestimated  the  radial 
transport  associated  with  the  increased  potential  drop,  or  that  the  outward  transport  occured 
over  a  limited  range  of  local  time  that  did  not  include  the  spacecraft. 


In  summary,  the  model  never  failed  to  predict  high  fluxes  when  they  were  observed. 
However,  it  sometimes  predicted  high  fluxes  when  they  were  not  observed,  and  it  failed  to 
predict  flux  dropouts.  On  one  occasion  it  predicted  a  dropout  during  a  low  flux  level.  The 
model  fluxes  for  prestorm  conditions  were  higher  than  reported  by  the  satellites. 
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Figure  4. 1  The  MSM  input:  Kp 


•  • 
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Figure  4.3  The  MSM  input:  Solar  wind  velocity  [kna/s] 
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Figure  4.6  The  MSM  input:  Polar  cap  potential  pattern 
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of  the  spatial  interpolation. 
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Figure  4.9  The  40keV  electron  differential  fluxes  of  the  MSM, 
Garrett  model  and  the  observations  by  the  satellite- 1 . 
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Figure  4.10  The  40keV  electron  differential  fluxes  of  the  MSM, 
Garrett  model  and  the  observations  by  the  sateUite-2. 

[  logio(  electrons  /cm^  /s  /sr  /keV  )  ] 
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Figure  4.1 1  The  40keV  electron  differential  fluxes  of  the  MSM, 
Garrett  model  and  the  observations  by  the  sateUite-3. 
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Figure  4.12  The  65keV  electron  differential  fluxes  of  the  MSM, 
Garrett  model  and  the  observations  by  the  satellite-1. 
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Figure  4.13  The  65keV  electron  differential  fluxes  of  the  MSM, 
Garrett  model  and  the  observation  by  the  satellite-2. 
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Figure  4.14  The  65keV  electron  differential  fluxes  of  the  MSM, 
Garrett  model  and  the  observations  by  the  satellite-3. 

[  logio(  electrons  /cm-  /s  /sr  /keV  )  ] 


6/30/90 


Section  5,  Page  1 


5.  MSM  Source  Code  Description  and  Documentation 


5.0  Language  and  Documentation 

The  MSM  code  was  developed  primarily  on  an  IBM  308  Id  and  is  in  ANSI  standard 
Fortran-77.  Before  delivery  to  the  Air  Force,  it  was  converted  for  operation  on  a  Microvax 
n  to  be  DEC  compatible.  The  B-field  matrices  were  computed  on  an  Apollo  10000. 

All  code  was  written  with  AFWL-TR-85-26,  Fortran-77  Computer  Program 
Structure  and  Internal  Documentation  Standards  for  Scientific  Applications,  as  a  guide. 
The  code  itself  is  internally  documented  and  has  been  carefully  maintained  throughout 
development. 

As  requested  by  the  Air  Force,  drafts  of  specific  sections  of  the  documentation 
consistent  with  DoD-STD-7935A  can  be  found  in  the  Appendix  of  this  report.  This 
documentation  was  placed  in  the  Appendix  because  the  paragraph  numbering  system 
described  by  DoD-S'ID-7935A  is  not  consistent  with  that  used  for  the  body  of  this  report. 
A  hard-copy  program  listing  and  tape  copy  are  included  with  this  report. 
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6.  The  Field  Line  Tracing  Program,  MAPINT 
6.0.  Background 

The  MSM  generates  values  for  the  electron  and  ion  fluxes  at  each  grid  point  on  a 
latitude  -  magnetic  local  time  grid  in  the  ionosphere  and  the  point  at  which  the  field  line 
from  each  grid  point  passes  through  the  equatorial  magnetic  reference  surface.  New  flux 
values  are  produced  every  15  minutes  of  magnetosphere  time.  MAPINT  is  the  program 
developed  by  Rice,  at  the  request  of  the  Air  Force,  to  meet  the  practical  requirement  for 
obtaining  the  particle  flux  at  any  arbitrary  satellite  position,  at  any  arbitrary  time.  MAPINT 
interpolates  and  extrapolates  the  MSM  output  to  relevant  times  and  locations.  This  is 
accomplished  by  tracing  magnetic  field  lines  and  using  the  interpolated  flux  from  the 
equatorial  plane  as  the  new  flux  value. 

MAPINT  is  interactive,  providing  operator  prompts  and  responding  with  the  results 
on  screen  displays.  MAPINT  can  also  form  the  basis  for  more  routine  or  sophisticated 
displays  of  MSM  output  such  as  those  discussed  in  Quarterly  Report  No.  1 1  and  presented 
by  Rice  at  the  Febmary  Review  meeting. 

MAPINT,  as  delivered,  is  a  core  program  that  uses  machine  specific  file  access  for 
it's  input  data  files.  The  delivered  version  is  self-contained  i.e.  has  it's  own  data  files,  and 
will  need  to  be  revised  for  use  with  the  MSM  on  any  new  machine. 

6.1  How  MAPINT  Works 

Since  the  MSM  magnetic  field  model  is  not  computed  "on-the-fly"  but  rather  uses  a 
table-lookup  of  precomputed  matrices,  the  tracing  of  the  magnetic  field  line  passing  through 
an  arbitrary  satellite  position  to  the  equator  cannot  be  done  directly.  Instead,  the  values  of 
the  B-field  input  parameters  (standoff,  Dst,and  lower  latitude  ^undary  of  the  aurora) 
nearest  in  time  to  the  desired  flux  time,  for  which  B-field  matrices  exist,  are  used  to 
generate  eight  field  line  traces  back  to  the  equator.  (The  delivered  version  of  MAPINT  is 
capable  of  including  both  dipole  tilt  and  collapse  configurations,  in  which  case  Liirty-two 
field  line  traces  are  performed.)  Figures  6. 1  and  6.2  illustrate  the  eight  nearest  neighbors  in 
parameter  space  and  the  tracing  of  the  resulting  eight  lines. 

Next,  the  satellite's  actual  field  line  crossing  point  is  determined  by  interpolating  the 
eight  model  crossing  points  using  weighting  factors  from  the  actual  parameter  space 
location.  Finally,  the  flux  at  the  satellite  point  is  taken  to  be  the  same  as  the  flux  at  the  point 
determined  above  as  interpolated  from  the  fluxes  at  the  adjacent  four  grid  points. 

6.2  MAPINT  Operating  Instructions 

The  following  is  a  list  of  program  input,  output,  program  checks  and  required  files 
for  MAPINT: 

Program  input: 

1 .  Time  (year,  Julian  day,  hour,  minute,  second)-integer  format. 

2.  Energy  (KeV)-real  format. 

3.  Particle  type  (electron,  H'*',  0‘*')-integer  format. (integer  code  included  in  prompt) 

4.  Point  in  space  from  which  to  start  mapping  e.g  a  satellite  position 

(GSM  coordinates,  Re)-real  format. 
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Program  output: 

1.  Spatial  spread  of  equatorial  reference  surface  points  used  for  interpolating  magnetic 
field  equatorial  crossing  point  (GSM  coordinates.  Re). 

2.  Point  on  the  magnetic  equatorial  plane  that  the  satellite  point  maps  to 
(GSM  coordinates.  Re). 

3.  Flux  at  the  input  data  point  (logbaselO(#/cm2-s-KeV-ster)). 

4.  Error  at  geostationary  orbit,  if  available. 

5.  Flag  telling  whether  fuU-traceback  or  KP-based  fluxes  were  used. 

Program  checks: 

1.  Input  time  is  reasonable  (ie,  hours<24,  minute<60). 

2.  Energy  >0. 

3.  Input  spatial  point  is  valid  (ie,  >1  Re). 

4.  Input  spatial  point  is  not  on  an  open  field  line. 

5.  Input  spatial  point  is  within  the  model  calculation  boundary. 

6.  Particle  type  is  correct. 

7.  Energy  and  particle  type  were  calculated  in  the  MSM  run. 

Files  needed  to  run  program  and  unit  numbers: 

1.  FLUX,  unit=15. 

2.  UPDAT,  unit=44. 

3.  VM,  unit=16. 

4.  BNDLOC,  unit=17. 

5.  XMIN,  unit=18. 

6.  YMIN,  units  19. 

7.  ZMIN,  unit*20. 

8.  The  magnetic  field  matrices  associated  with  the  run. 

The  program  is  fully  interactive.  It  requests  the  input  information  and  displays  the  results 
on  the  terminal.  The  following  pages  are  copies  of  an  actual  session: 


DnSLIO710l  Execution  begins... 

ENTER  LAST  2  DIGITS  OF  VERB  AND  JULIAN  DRV 

7 

86  112 

ENTER  HOUR,  NINUTE,  AND  SECONDS 

7 

18  10  30 

INPUT  ENERGY  (KEU)  AND  1  FOR  ELECTRONS,  2  FOR  H+, 
OR  3  FOR  0+ 

7 

11.5  1 

ENTER  SPATIAL  POINT  OF  INTEREST  IN  GSN  (RE):  X,V,Z 

7 

0  -6.5  1.15 

HAPPING  POINTS  USED  FOR  INTERPOLATION  UERE  SPREAD: 
0.273001286E-01  RE  ALONG  THE  XGSH  DIRECTION 
0.657811513E-01  RE  ALONG  THE  VGSH  DIRECTION 
0.261033997E-01  RE  ALONG  THE  ZGSH  DIRECTION 

INTERPOLATED  HAPPING  POINT:  XGSH  VGSH  ZGSH 
-0,229525117E-01  -6.85171817  -0 .299130827E-01 


LOG  FLUX  AT  THE  SATELLITE  =  1.12708187  HOLDING  RICECSUH 


TRACE  FROH  ANOTHER  STARTING  POINT 

UlTH  SRHE  TIHE  PRRRHETERS?  (1  =  VES) 

,  7 

:  1 

;  INPUT  ENERGY  (KEU)  AND  1  FOR  ELECTRONS,  2  FOR  H+, 
i  OR  3  FOR  0+ 

7 

:i50  I 
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,  7 

!5. 16  3.96  1  .15 

iHAPPING  POINTS  USED  FOR  INTERPOLATION  UERE  SPREAD: 
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'  0.152785192E-01  RE  ALONG  THE  VGSH  DIRECTION 
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/ 
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0 
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7 

1 

I  ENTER  LAST  2  DIGITS  OF  VEAR  AND  JULIAN  DRV 
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7 

88  112 
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7 

18  23  15 
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I  ’ 
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2 
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7 
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6.3  MAPINT  Validation  Using  the  April  21-23,  1988  Event 

As  a  validation  and  demonstration  exercise,  MAPINT  was  applied  to  a  hypothetical 
spacecraft  located  at  3  Re  and  45  degrees  inclination  (GSM  coordinates:  X=~2.12,  Y=0, 
Z=2.12)  for  three  times  during  the  April  21-23,  1988  event.  The  times  selected  were 
intended  to  be  illustrative  of  expected  behavior  of  the  magnetic  field  during  certain  phases 
of  the  storm.  The  times  were  18:30  on  day  1 12,  and  00:00  and  06:00  on  day  1 13.  Figure 
6.3  shows  the  results. 

In  figure  6.3,  we  have  a  cartoon  of  the  field  line  shapes  and  positions  for  each  of  the 
three  times,  for  field  lines  originating  at  the  spacecraft.  The  field  lines  are  not  drawn  to 
scale.  The  Log  of  the  electron  fluxes  and  the  equatorial  crossing  points  of  the  satellites,  as 
computed  by  MAPINT,  are  shown.  We  also  show  the  position,  in  the  storm,  of  each  time 
as  seen  in  the  geophysical  parameter  diagram  at  the  bottom  of  the  figure. 

The  first  time,  18:30  on  day  112,  is  prestorm  -  early  growth  phase.  The  figure 
shows  that  the  field  is  reasonably  dipolar,  in  other  words  not  stretched  tailward  too  much 
by  tail  currents.  The  second  time,  the  start  of  day  113,  shows  a  highly  stretched  field  line 
expected  of  the  expansion  phase,  before  the  collapse  or  injection  has  reached  the  L-shell  of 
the  spacecraft.  The  last  time,  06:00  on  day  113,  shows  the  field  line  connecting  the 
spacecraft  to  have  coUapsed  to  the  more  dipole-like  configuration,  as  would  be  expect^  for 
late  injection. 

The  magnetic  field  morphology  and  dynamics  exhibited  here  are  generally  what  are 
expected.  Since  there  are  no  measurements  of  the  magnetic  field  available,  there  can  be  no 
true  validation  beyond  this  qualitative  analysis.  Moreover,  we  have  no  satellite  particle  flux 
data  at  the  hypothetical  location.  Validation  of  the  B-field  model  is  actually  accomplished 
by  validation  of  the  MSM,  as  discussed  in  sections  3  and  4.  We  are  satisfied  that  MAPINT 
is  providing  a  correct  extrapolation  of  the  MSM  flux  values  and  is  using  correctly  the  B- 
field  model  available  to  it. 

For  acceptance  purposes,  entry  of  the  given  times  and  satellite  location  in  MAPINT 
should  reproduce  the  fluxes  and  field  line  equatorial  crossing  points  shown  in  figure  6.3. 

Detailed  documentation  on  MAPINT  can  be  found  in  the  Appendix. 
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6.3  MAPINT  Validation  Using  the  April  21-23,  1988  Event 

As  a  validation  and  demonstration  exercise,  MAPINT  was  applied  to  a  hypothetical 
spacecraft  located  at  3  Re  and  45  degrees  inclination  (GSM  coordinates:  X=-2.12,  Y=0, 
Z=2.12)  for  three  times  during  the  April  21-23,  1988  event.  The  times  selected  were 
intended  to  be  illustrative  of  expected  behavior  of  the  magnetic  field  during  certain  phases 
of  the  storm.  The  times  were  18:30  on  day  1 12,  and  (X):00  and  06:(X)  on  day  113.  Figure 
6.3  shows  the  results. 

In  figure  6.3,  we  have  a  cartoon  of  the  field  line  shapes  and  positions  for  each  of  the 
three  times,  for  field  lines  originating  at  the  spacecraft.  The  field  lines  are  not  drawn  to 
scale.  The  Log  of  the  electron  fluxes  and  the  equatorial  crossing  points  of  the  satellites,  as 
computed  by  MAPINT,  are  shown.  We  also  show  the  position,  in  the  storm,  of  each  time 
as  seen  in  the  geophysical  parameter  diagram  at  the  bottom  of  the  figure. 

The  first  time,  18:30  on  day  112,  is  prestorm  -  early  growth  phase.  The  figure 
shows  that  the  field  is  reasonably  dipolar,  in  other  words  not  stretched  tailward  too  much 
by  tail  currents.  The  second  time,  the  start  of  day  113,  shows  a  highly  stretched  field  line 
expected  of  the  expansion  phase,  before  the  collapse  or  injection  has  reached  the  L-shell  of 
the  spacecraft.  The  last  time,  06:00  on  day  113,  shows  the  field  line  connecting  the 
spacecraft  to  have  collapsed  to  the  more  diople-like  configuration,  as  would  be  expected  for 
late  injection. 

The  magnetic  field  morphology  and  dynamics  exhibited  here  are  generally  what  are 
expected.  Since  there  are  no  measurements  of  the  magnetic  field  available,  there  can  be  no 
true  validation  beyond  this  qualitative  analysis.  Moreover,  we  have  no  satellite  particle  flux 
data  at  the  hypothetical  location.  Validation  of  the  B-field  model  is  actually  accomplished 
by  validation  of  the  MSM,  as  discussed  in  sections  3  and  4.  We  are  satisfi^  that  MAPINT 
is  providing  a  correct  extrapolation  of  the  MSM  flux  values  and  is  using  correctly  the  B- 
field  model  available  to  it. 

For  acceptance  purposes,  entry  of  the  given  times  and  satellite  location  in  MAPINT 
should  reproduce  the  fluxes  and  field  line  equatorial  crossing  points  shown  in  figure  6.3. 

Detailed  documentation  on  MAPINT  in  accordance  with  DoD-STD-7935A  is  given 
in  the  Appendix. 
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7.  Summary 

Rice  has  successfully  developed  a  practical  computer  model  capable  of  specifying 
electron  and  ion  fluxes  in  the  magnetosphere  during  geomagnetic  storms,  based  on  input 
data  from  the  Environmental  Database. 

The  Magnetospheric  Specification  Model  represents  the  first  effort  at  a  large-scale 
computational  model  of  the  Earth’s  magnetosphere  that  is  designed  for  operational  use. 
The  MSM  computer  algorithm  represents  a  substantial  effort  to  bring  the  observational 
knowledge,  theoretical  understanding,  and  computer-modeling  technology  that  has  been 
developed  in  the  last  thirty  years  of  space  research  to  bear  on  practical  operational 
problems. 

Our  testing  of  the  MSM,  as  presented  in  Sections  3  and  4  of  this  Repon,  indicates 
that  it  has  achieved  significant  success  in  its  primary  goal  of  representing  fluxes  of  kilolvolt 
electrons  in  the  geostationary-orbit  region.  It  clearly  provides  good  and  useful 
specifications  of  these  fluxes,  and  it  represents  a  clear  advance  in  the  state  of  the  art. 

The  application  program  MAPINT,  which  is  provided  with  the  MSM,  provides  the 
capability  for  mapping  fluxes  from  an  arbitr^  point  in  the  middle  magnetosphere  to  the 
equatorial  plane,  thus  providing  the  capability  for  flux  specification  for  non-equatorial 
spacecraft. 

A  major  advantage  for  the  future  is  that  the  MSM  calculates  particle  fluxes  within  the 
framework  of  a  general  magnetospheric  model  that  consistently  calculates  most  other  large- 
scale  parameters  of  the  physical  system,  including  the  fluxes  of  precipitating  auroral 
electrons  and  ionospheric  electric  fields.  This  gives  the  MSM  great  potential  for  growth  in 
terms  of  wider  and  more  valuable  specifications. 

A  major  limitation  of  the  present  MSM  is  that  many  of  its  capabilities  for  calculating  a 
wide  variety  of  magnetospheric  parameters  could  not  be  tested  as  extensively  as  we  would 
have  liked  against  real-time  observational  data.  This  limitation,  which  was  due  to  the 
configuration  of  available  spacecraft  and  instrumentation,  should  disappear  in  large  pan 
over  the  next  few  years,  as  new  observational  data  become  available  to  us.  The  testing  of 
the  MSM  that  is  planned  for  the  follow-on  contract  will  provide  much  better  quantitative 
information  on  the  model's  capabilities  for  specifying  parameters  beyond  the  basic  35-100 
keV  geosynchronous  electrons.  Additional  particle  detectors  that  are  now  being  installed 
on  DoD  geostationary  spacecraft  will  provide  checks  on  predictions  of  low-energy 
particles.  The  CRRES-SPACERAD  spacecraft  will  provide  a  rich  and  invaluable  source  of 
data  for  L  <  6.6.  Data  that  are  now  being  synthesized  from  the  great  magnetic  storm  of 
March  1989  will  allow  us  to  verify  and  quantify  the  MSM’s  accuracy  for  the  most  extreme 
geomagnetic  conditions.  The  flexible,  modular  structure  of  the  MSM  will  allow  us  to 
carry  out  simple  "fixes"  by  changing  very  few  lines  of  code,  if  the  extended  testing 
indicates  the  need  for  such  changes. 

In  summary,  the  objectives  of  the  Magnetospheric  Specification  Model  have  been 
met,  and  a  model  algorithm  is  ready  for  adaptation  for  use  in  an  operational  setting,  where 
the  goal  is  real-time  and  retrospective  specification  of  hazardous  charged-particle  fluxes. 
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